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NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 


TECHNICAL NOTE D-1455 


THE N-BODY CODE - A GENERAL FORTRAN CODE FOR 

SOLUTION OF PROBLEMS IN SPACE MECHANICS 

BY NUMERICAL METHODS 

By William C. Strack, Wilbur F. Dobson, 
and Vearl N. Huff 


SUMMARY 

A general astronomical integration code designed for a large class of prob- 
lems in space mechanics that may be solved by numerical integration is described. 
The equations of motion provide for the effects of up to eight gravitating celes- 
tial bodies, oblateness and aerodynamic forces from the celestial body at the 
problem origin, propulsion system thrust, and rotation of the body at the origin. 


INTRODUCTION 

The general problems of space mechanics (i.e. , n-bodies plus nonconservative 
forces such as thrust) cannot be solved analytically. Therefore, numerical inte- 
gration through the use of computing machinery is usually employed. 

Several codes have been written for the numerical solution of problems in 
orbit mechanics; for example, the Themis Code of reference 1 is a double- 
precision code intended primarily for close satellites or interplanetary coasting 
flight. Reference 2 describes a space-trajectory program of considerable merit. 

A listing of several other trajectory codes may be found in reference 3. 

The gener a] purpose code described herein has several distinctive features 
not all of which are found in any one of the previously available codes. As de- 
scribed herein, this code is designed to operate on an IBM 704 computer that has 
an 8000 word (8 K) memory and at least 1 K of drum. The fact that the program is 
written in FORTRAN should make it applicable to installations having other types 
of equipment that accept the FORTRAN language. An edition of this program (dif- 
fering primarily in that segmenting of the program is not required) is available 
for an TPM 7090 computer that has a 32-K core. 

The program is compartmented into 25 subroutines to facilitate modifications 
for specific problems. The integration is carried out in either rectangular co- 
ordinates or orbit elements at the option of the user. A compact ephemeris that 



occupies about one-seventh of a reel of tape is utilized for positions and veloc- 
ities of the planets (except Mercury) and the moon. An atmosphere is included so 
that aerodynamic forces may be considered. 

STATEMENT OF PROBLEM 

The problem to be solved may be stated as follows: Given certain initial 

conditions , compute, using three degrees of freedom, the path of an object, such 
as a space vehicle, subject to any or all of the following forces: 

origin body gravitational field 

other celestial body gravitational fields 

propulsive thrust 

aerodynamic forces 

any other defined forces 

or, in equation form, with respect to a noninertial Cartesian coordinate system, 



where n equals the number of perturbating bodies and 7 denotes the del oper- 
ator. (All symbols are defined in appendix A.) 


Origin Body Gravitational Field and Oblateness Perturbations 


The first term, VQ, in the equation of motion (eq. (l)) represents the 
gravitational forces due to the origin body. When the origin body is spherical 
and made up of homogeneous layers, this term becomes simply -pr/r^. In the case 
of the Earth, however, the effect of oblateness may be important, and additional 
terms must be added to account for the oblateness effects. The expression for 
the gravitational potential U of an oblate spheroid may be written, according 
to reference 4, as 



where the x,y plane lies in the equatorial plane. The components of gravita- 
tional acceleration are as follows : 
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( 3 ) 


The first terms exist for a spherical planet composed of concentric layers of 
uniform density. The terms containing J and K are frequently called the sec- 
ond and fourth harmonic terms, while J and K are known as the harmonic coef- 
ficients. 

It is expected that ohlateness perturbations need to he computed only for 
the origin tody, since at large distances, such as that between celestial bodies, 
the gravitational field of an oblate body is essentially an inverse-square field. 
Consideration of oblate bodies other than the Earth requires only different val- 
ues of J and K if that body's rotation axis is parallel to the z-axis. When 
the body has triaxial asymmetry or when the z-axis cannot conveniently be alined 
with the rotation axis of the origin body, the equations must be extended if ob- 
lateness is to be included. 


Celestial Body Perturbations 

The presence of more than one gravitating body in addition to the object re- 
sults in the inclusion of the second term of equation (l). The evaluation of 
this term requires a knowledge of the positions of the bodies as a function of 
time. The degree of precision desired determines the method to be used to obtain 
the positions such as elements of ellipses or an ephemeris. 


Propulsive Thrust 

The propulsive acceleration is completely specified by a direction and a 
magnitude. The thrust direction may be referred to the velocity vector by two 
angles: a, the angle between the velocity and the thrust vectors, and |3, the 
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angle between the orbit plane and the velocity-thrust plane. The sense of each 
angle is indicated in sketch (a). 



y 


(a) 

The velocity may be referenced with respect to one of several coordinate 
systems. If the computation refers to a takeoff of a rocket or winged vehicle, 
the coordinate system rotating with the Earth may be preferred. In such cases 
the relative velocity (i.e., the velocity of the object relative to the atmos- 
phere) will serve to orient the thrust vector. Resolution of the thrust-vector 
components along the x,y, z axes is shown in appendix B. 

The thrust magnitude of a rocket engine is 

F - mlg c - PA e (4) 

This relation is sufficient for many space powerplants and is used in the present 
program. 


Aerodynamic Forces 

The aerodynamic forces are usually divided into the two components, lift and 
drag. The drag force is directed opposite to the relative wind vector, and the 
lift vector is perpendicular to the relative wind vector. The angles a and 0, 
defined in the previous section, serve as the angles of attack and roll, respec- 
tively. Yaw effects are not considered. Resolution of the lift and drag vectors 
into components along the x,y,z axes is given in appendix B. 

The magnitudes of the lift and drag forces may be conveniently determined 
through use of a tabular group of coefficients in relatively simple equations. 

The lift and drag magnitudes may then be expressed (as is usual in aerodynamics) 
as 
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( 5 ) 

( 6 ) 


L - C^(a,l^)qS 
D = Cp(Cj* ji ,N^)qS 


where 

a = a(t) 

C T = f- (K-)sin a 
L I'M 

C D = C D,0 + C D,i = C D,0^^ + 
q =|p(V') 2 
p = p(P,T) = p(h) 

If a(t), Cp q(%)> and f 2(%) axe assume( i to be quadratic functions and 

p is assumed to be constant, the expressions for a, P, Cp q and Cp ^ 
become 


a 


a ll + a 12 t + a l ,3 t2 


P - P 0 

C L = ( a 21 + a 22 W M + a 23 1 ^d) S:Ln a 
C D,0 = a 31 + a 32 W M + a 33 H M 
C D,i = ( &41 + a 42% + & 43 N m) C L 

where the quadratic constants j may have different values for different re- 
gions of the independent variables t and Hy. 


It should be remembered that these choices are arbitrary and are not re- 
strictive because other functions may easily be used by simply changing the equa- 
tion where it appears in the program. In fact, any propulsion system and aerody- 
namic configuration can presumably be incorporated by writing proper thrust and 
aerodynamic subroutines . 

The pressure, temperature, and density may be determined as a function of 
altitude in accordance with the ICAO standard atmosphere. The oblate Earth model 
is used to determine the altitude. 
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Other Forces 


The X forces may he any forces such as electrostatic, magnetic, or solar 
radiation pressure that affect the trajectory. 'While these forces are not con- 
sidered further herein,, their inclusion would usually he feasible and would he 
similar to thrust, lift, and drag. 


METHOD OF SOLUTION 

The method of solution selected for the stated problem is presented in this 
section. A later section discusses the FORTRAN coding. 

A description of several numerical integration techniques and their relative 
merits are contained in reference 5, A straightforward method for finding the 
position of the object as a function of time is to integrate the total accelera- 
tion of the object expressed in rectangular components. An example of this 
method is Cowell T s method (ref. 5). 

However, when the system under investigation consists of two nonoblate 
bodies (one of which is the object) with no forces other than gravitational at- 
traction forces, an exact analytical solution for the motion of the body exists. 
Further, if the conditions of the actual problem are such as to approximate the 
two-body problem closely, another approach is to use the exact two-body solution 
as a basis and simply integrate the changes in the two-body parameters, since 
they should be slowly varying. This technique, sometimes called the "variation 
of parameters," will be referred to as "integration of orbit elements." 

Since problems both remote and near to the exact two-body problem are en- 
countered in orbit mechanics, and since either type of problem is solved more ef- 
ficiently by using the technique most suitably applicable, it was considered de- 
sirable to use either of the previously mentioned integration techniques at will. 
Accordingly, two methods of integration are provided in the program, namely, rec- 
tangular coordinates and orbit elements. 


Integration Variables 

In order to use either of these integration techniques, it is necessary to 
select a suitable set of variables for integration. Because a differential equa- 
tion may determine the mass of the object (i.e,, spacecraft), mass has been se- 
lected as a variable to be integrated. Selection of the remaining parameters 
follows in the subsequent paragraphs. 

Rectangular coordinates . - In the first technique, the total acceleration 
components X,y, and V are integrated to obtain x,y, and z where x,y, and z are 
the rectangular components of the origin-to-object radius r. The positive 
x-axis points in the direction of the mean vernal equinox of 1950.0. The posi- 
tive y-axis lies in the mean equator of 1950.0 and is perpendicular to and coun- 
terclockwise from the positive x-axis. The z-axis points north and completes the 
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righthanded orthogonal set. The integration in rectangular coordinates involves 
numerical solution of three second-order linear differential equations; that is, 
a double integration is required for integrating the accelerations to obtain ve- 
locities and the velocities to obtain positions. The rectangular variables have 
advantages of complete generality and a minimum amount of computing per step. 

Orbit elements . - In the variation-of -parameters technique, a set of six 
independent two-body parameters called orbit elements are integrated. These six 
parameters may be arbitrarily chosen from a host of possibilities- The set se- 
lected for this program is composed of the eccentricity e, the argument of peri- 
center cd, the equatorial longitude of ascending node ft, the inclination of the 
orbit plane to the equatorial plane i, the mean anomaly M, and the semilatus 
rectum p. The transformation equations between the two sets of variables are 
given in appendix C. 

The integration of orbit elements requires the numerical solution of six 
first-order linear differential equations. The rather involved transformation by 
which the three second-order linear differential equations in X,y,^ are reduced 
to six first-order equations in 6, d>, ft, i, M, and p is contained in refer- 
ence 6. Integration in orbit elements is frequently advantageous because the 
smaller orbit -element derivatives may permit larger integration intervals that 
result in fewer steps. In the special case of two-body motion, the derivatives 
are zero (except M, which is a constant). 

Mathematical difficulties may arise occasionally with most sets of orbit 
elements. In particular, for the selected set, these occur when e approaches 
unity (parabolic trajectory), which causes a loss of numerical accuracy in the 
frequently used quantity (l - e 2 ); and when an asymptote is approached too 
closely, which causes numerical difficulties in the iterative solution for eccen- 
tric anomaly from Kepler* s equation. The selected solution to these difficulties 
is to shift temporarily to rectangular-coordinate integration whenever the diffi- 
culty arises. 


Integration Method 

It Is clear that regardless of the choice of integration technique, the mag- 
nitudes of the derivatives of the variables to be integrated may vary consider- 
ably along the trajectory. With fixed step size (constant intervals in time), 
the integration scheme will take unnecessary steps in the regions where the 
changes in the derivatives are small and thus will waste computing time and In- 
crease roundoff error. When the derivatives are large and change rapidly, a 
fixed step size will result In large truncation error (error due to excessive 
step size). Thus, in the interest of computing accuracy and economy, use of 
variable step size along the trajectory becomes desirable. 

One of the integration schemes that allows variable step-size control to be 
incorporated easily is the Runge-Kutta scheme. For this and other reasons, it 
was decided to use a fourth-order Runge-Kutta method with variable step-size 
control. 
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Truncation error and step size may be controlled by examining the relative 
errors betveen the fourth-order Runge-Kutta integration scheme and a lover -order 
integration procedure. The arbitrarily chosen lov-order integration scheme vas 
an unequal-interval Simpson rule method. Details of the fourth-order Runge-Kutta 
integration method and the step-size control are given in appendix D. Roundoff 
error may be reduced by accumulating the integration variables in double preci- 
sion. 


Origin Translation 

As noted previously, machine computing time and roundoff error may be mini- 
mized by maximizing the integration interval. The largest intervals are possible 
in orbit elements vhen the celestial body at the problem origin is the one that 
has the greatest influence on the vehicle motion. For this and sometimes other 
reasons, it may become desirable to translate the problem origin occasionally as 
the vehicle moves along its path. 

Such translations of the origin may be made vhen the object enters a body's 
"sphere of influence," that is, the sphere about a body vithin vhich the greatest 
influence upon the object is due to forces originating from that particular body. 
In this program, the orientation of the coordinate system is alvays alined vith 
the system determined by the Earth's mean equator and equinox of 1950.0, as is 
standard in astronomy. 


THE CODE AND ITS USAGE 

The stated problem vas programmed in FORTRAN routines that are separately 
designed to accomplish one task but vhen combined form a complete program. This 
feature facilitates modifications. 

The program is labeled as a general-purpose code, but an efficient general- 
purpose code cannot be a reality. As a result, this code is not especially gen- 
eral, but an attempt has been made to retain efficiency and to provide for easy 
modification of the routines to recover generality as needed. For example, the 
program is an "open system," that is, it solves an initial value problem. There 
is no link provided to obtain specific end conditions. Provision of this link is 
left to the user for his specific needs. In particular, vhen certain end condi- 
tions of a trajectory are to be met by determining the correct initial conditions 
(tvo-point boundary value problem), the user may program an iteration scheme to 
compute initial conditions from end conditions of previous runs. 

The code is designed to operate on an IBM 704 computer that has an 8-K core 
and drum and also a number of tape units. To operate the code on an 8-K com- 
puter, it is necessary to divide the program into tvo segments (core loads). The 
program of segment 1 arranges certain data in the core. The program of segment 2 
overwrites the program but not the data of segment 1 vhen it is called for. Fig- 
ure 1 is a simplified diagram that shovs hov the various major subprograms are 
arranged in the segments. The segmenting vas done as efficiently as possible in 
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Figure 1. - Block diagram of program segments shoving principal subprograms . 







terms of execution time, but further gains can be realized by users of larger 
computers vho may wish to modify the code to utilize the increased computer ca- 
pacity. 

In the following sections, discussing the program in terms of the FORTRAN 
variables and routines is sometimes desirable. A glossary of these variables 
is given in appendix E. 


Ephemerides 

To determine the position of each celestial body, there is offered a choice 
between ellipses and a precision ephemeris. Any appropriate ellipse data may be 
used, and an example of such data is given in table I. 

The precision-ephemeris tape that is used in the program was so made that 
position and velocity were obtainable through the use of a fifth-order polynomial 
whose coefficients are stored on tape. The details concerning the making of the 
tape and its structure are given in appendix F. This master tape Is a merged 
ephemeris containing all the planets (except Mercury), the moon, and the Earth- 
moon barycenter from October 25, 1960 to about 2000 (except for the moon, which 
has an ending date of 1970). The Earth ephemeris is called ,T sun ir because it 
gives sun to Earth distances. 

Direct use of the master merged ephemeris tape would, in general, be waste- 
ful of computing time, since excess tape handling would occur in order to bypass 
data not required for the particular problem. To minimize tape handling during 
execution, a shorter merged ephemeris containing only that data needed for a spe- 
cific problem is constructed at execution time. Several of these working ephem- 
erides may be constructed before the integration of the problem, (Several prob- 
lems may be loaded simultaneously with the same ephemeris, or each problem may 
require a distinct ephemeris, or several ephemerides may be desired for a single 
problem. ) 


Step Size and Output Control 

Truncation error and step size are controlled by computing the relative 
errors between the Runge-Kutta integration and the lower-order integration proce- 
dure. If the greatest relative error between the methods is greater than a maxi- 
mum limit (ERLIMT), the integration step will be repeated after a smaller step 
size is computed. In either case, a new step size is computed from the relative 
errors of the previous steps and is intended to result in an error that is close 
to a reference value (EREF). Further, the step size may then be reduced by the 
output controls. In any case, a step can be no larger than three times the size 
of the previous successful step. (See appendix D. ) 

Output is sometimes desired at specific points along the trajectory, while 
at other times this is unimportant. This option is provided for the user so that 
he may choose output to occur at equal intervals in step number or equal time 
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intervals (which places a constraint on the step size). Also, he may choose to 
change from one mode to another along the trajectory. These choices of output 
spacing are effected through the use of the FORTRAN variables MODOUT, DELMAX, 
STEPS, and TMIN, which is explained under the MODOUT entry of table II, a table 
of program control parameters. 


Computer Output 

A basic output format was programmed to serve as a basis for modification 
and is illustrated in table III. It is Intended that a user of the code modify 
the output to suit his purpose. In addition to examining the normal output, it 
Is sometimes desirable to examine the error-control data, such as the relative 
errors in the integration variables, along the path. These data are printed as a 
single block after completion of the problem if the sign of the input error ref- 
erence value EREF is negative. The sign of EREF is Irrelevant in the error- 
control portion of the program since its absolute value is taken. 


Computer Input 

The user has a choice of three possible sets of Input data that specify po- 
sition and velocity: (1) the six orbital elements, (2) the three Cartesian com- 

ponents of both velocity and position, and (3) the latitude, longitude, azimuth, 
elevation, velocity, altitude, and time. 

The third set mentioned is programmed for the Earth only where the latitude 
and longitude are the geocentric latitude and longitude measured from the equator 
and Greenwich, respectively. The azimuth angle is measured in a plane tangent to 
the sphere of radius r at the point on the sphere determined by the geocentric 
latitude and longitude, and relative to the local meridian, positive eastward 
from north. The elevation angle is then measured in a plane normal to the tan- 
gent plane, positive outward (sketch (b)). The tangent plane is taken horizontal 
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Long - longitude measured from 
Greenwich in earth’s 
equatorial plane, posi- 
tive east 

Lat » latitude, measured posi- 
tive north, geocentric 

Azi = azimuth angle, measured 
east from north from 
local meridian 

Elev - elevation angle, positive 
outward 

Vel « vehicle’s initial velocity 

r - radius of vehicle from 

Earth’s center 


n 



with the effects of oblateness and rotation considered if these effects are "on. " 
If oblateness and rotation are "off," the horizontal is perpendicular to the ra- 
dial direction. This input option ignores the correction between universal time 
and ephemeris time and between the instantaneous equator and equinox and those cf 
1950.0. 

A list of input instructions is contained in appendix G along with an input 
check list. 

The input routine described in reference 7 was used because of its simplic- 
ity; however, another input routine may be used if it is desired. 


Sequence of Operations 

Before the program begins to integrate a trajectory, it performs an assort- 
ment of operations that may be called "initialization. " All these operations ar e 
expected to be done once or only a few times during the trajectory integration 
and, for this reason, are contained wholly in segment 1. Likewise, at the end of 
a problem, a return to the segment 1 causes several concluding operations to be 
performed. A condensed description of the operations carried out In segment 1 is 
contained in the flow diagram of figure 2. Other than the normal end of a prob- 
lem (reaching a maximum number of integration steps or a particular time) there 
is only one way in which segment 1 may be called by segment 2, namely, a trans- 
lation of the origin. When the translation occurs, segment 1 is needed to re- 
order the body list and perhaps to cause input or ephemeris change. 

After completion of the initialization, which leaves numerical data stored 
in the common area, segment 1 is overwritten by segment 2, which may be termed 
the integration segment. 


CODING 

General 

Appendix H contains the code listing of the program. Although most of the 
program is coded in basic FORTRAN II, on several occasions it was preferable to 
use the pseudo-SAP statements of FORTRAN II. Typically, the pseudo-SAP statement 
LXD (l),I is used whenever the index I was to be transferred from one subroutine 
to another (since FORTRAN II does not do this automatically). Wherever such a 
statement appears, the FORTRAN II statement I = I can be used instead to accom- 
plish this initialization but with additional commands. 

Some of the FORMAT statements are of the G-type. These statements will 
print output in I, E, or F format depending on the nature of the variable. 
Fixed-point variables will take the I format, while floating-point variables will 
assume the F format unless the magnitude of the variable falls outside the useful 
F range, in which case the E format is used. FORTRAN facilities that do not ac- 
cept the G-type format statements may easily substitute E-type formats. 
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Figure 2. - Flow diagraa of segment 1. 






Table IV Is a map of COMMON allocation (blanks are left for the user) and 
table II contains a description of the program control parameters. The elements 
of the Integration variable array (XERIM) are given in table V. The assumed 
values of the astronomical constants are given In table VI. These values are 
easy to change to any set desired. A selected set is given in reference 8. 


Examples 

Two examples of code usage are presented in the following sections. The 
first example is a problem of raising a low-altitude satellite into a 24-hour 
orbit by using low-tangential acceleration. The other example is a more complex 
problem involving a ground-launched lunar probe with a three-stage rocket. Both 
problems were selected to illustrate the usage of the program rather than to at- 
tempt a detailed analysis of the example problem. 

Example I: Low-tangential thrust . - The trajectory to be determined is that 

used to raise a 3850-kilogram package from an initial 300-statute-mile circular 
equatorial orbit to a 24-hour orbit using a 60, 000-watt nuclear electric system 
with a specific impulse of 2540 seconds and an overall efficiency of 40 percent. 
The required engine parameters may be calculated as follows: 

thrust force: 


F = 


2 V) 

ig c 


2 X 60,000 X 0.4 
2540 X 9.80665 


= 1.927 newtons 


initial acceleration: 


_F_ 

“0 


1.927 

3850 


5 . 0051948 X 10" 4 m/sec 2 


propellant flow rate: 

-* - - 2 S 4 oWe 0665 = 7 - 7361935 X 10- 5 kg/sec 

A detailed account Is given in the following paragraphs for the solution of 
this problem by the prescribed program. Only those features of the program that 
have a direct bearing on this particular problem are discussed. Additional pro- 
gram features are discussed in the account of the second example problem. It may 
prove beneficial to refer to figure 2 during these two discussions. 

It is assumed in the program that all memory data stores are cleared (set 
equal to zero) before operation begins. Control begins when the routine MAIN 1 
is entered in segment 1. After several noninfluencing commands, the reading of 
a "clock" takes place at statement 10 and this value is stored. This value is 
later subtracted from the subsequent reading in order to yield the computing 
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time. (All references to the "clock" may be deleted without ill effect.) Then 
a set of so-called "standard data" is initialized by executing subroutine STDATA. 
Before initializing, STDATA clears most of COMMON C. 

The next step is calling for input at statement 21. The following list of 
parameters constitutes the input: 


Parameter 

FORTRAN 

name 

Value 

Initial mass, m q, kg 

RMASS 

3850 

Semilatus rectum, p, m 

P 

6.86X10 6 

Specific impulse, I, sec 

SIMP 

2540 

Flow rate, -m, kg/sec 

FLOW 

7. 7361935X10 -5 

Time limit, sec 

TMAX 

a 42605 

Initial step size, sec j 

DELT 

a 1500 

Step number limit, steps 

STEFMX 

a 2000 

Frequency of output, 
steps /output 

STEPS 

a 200 


a Assumed value. 


Variables such as eccentricity and mean anomaly that are initially zero are 
not included in this list since all memory data stores are initially zero. 

In accordance with the input routine of reference 7, the input cards may 
appear as 

$DAIA = 1 * STABLE ♦ 4 1 = RM ASS * 47 = P * 5 = S IMP * 33= $$ IDENTIFICATION AND 

FLOW*10 = DELT *30=TMAX * 20 = S TEPMX * 2 1 = STEP S/ TABLE DEFINITION 

RMASS=3850#SIMP=25A0*FLOW=7.7361935E-5 $$ VEHICLE MASS* ISP* MASS FLOW 
P = 6,86E6*TMAX = A2605 *STEPMX = 2000 $$ SEM I LATUS -RECTUM > TIME LIMIT* STEP LIMIT 

DELT=1500*STEPS=200 5$ INITIAL STEP SIZE* OUTPUT EVERY 200TH STEP 

SD Aj A= 1 » $$ LAST CARD 

where the entries between the $TABLE and slash (/) reference the subsequent en- 
tries to the second argument C of the calling statement. Thus, for example, 

RMASS is equivalent to C(4l), the 41 s ^ location from the beginning of COMMON C. 

Several commands follow the input none of which has an important effect on 
this particular problem with one exception; subroutine ORDER (part 11) computes 
the gravitational constants \x and n/jl. The initialization process is now 
completed. 

Segment 2 overwrites segment 1, except COMMON C(l) to COMMON C(800), and 
control begins when the routine MAIN 2 is entered. Immediately, the tape that 
stores the two segments (tape 2 at Lewis) is rewound to position this tape at the 
beginning of segment 1. 
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The next sequence is that of integrating the first two steps. These two 
steps are of equal size and are integrated before an error check is made. If the 
first two steps are satisfactory (determined by statement 25 ), the remaining 
steps are integrated while the relative error is being checked at the end of each 
step. Parts 1 and 5 of MAIN 2 are concerned solely with this starting phase. 

Part 1 sets up the starting sequence and causes the initial conditions to appear 
on the output sheet. Parts 2 to 4 accomplish the Runge-Kutta Integration for a 
single step. 

The derivatives used in the integration are obtained from subroutine EQUATE. 
The first half of this subroutine finds the Cartesian coordinates and velocities 
through use of Kepler *s equation. The thrust Is computed in statement 34, and 
then subroutine THRUST is called to determine the components of the thrust accel- 
eration in the Cartesian coordinate system. (After control is returned to sub- 
routine EQUATE, the thrust acceleration is resolved into circumferential, radial, 
and normal components. ) Finally the derivatives of the orbit elements are calcu- 
lated, and a return is made to MAIN 2. 

After the Runge-Kutta integration is performed, the error check is made in 
part 5B (part 6 after the starting sequence) by computing the difference between 
the Runge-Kutta integration and the low-order integration. Subroutine ERRORZ is 
called to determine the largest of the relative errors. If the largest of the 
relative errors is greater than the limit value, ERLJMT (set in STRATA), part 8, 
which computes a smaller step size for the same interval, is entered and control 
is returned to part 1. If the greatest relative error is smaller than the limit 
value, part 7, which advances the variables of integration, is entered and calls 
subroutine STEP to compute the next step size and print out the variables of the 
first step. Part 7 also counts the revolutions past the x-axis and adjusts the 
argument of pericenter and mean anomaly to within ±jt to retain accuracy in the 
sine-cosine routines. If the step size exceeds l/2 revolution, the revolution 
count may be short by an integral number. Control is finally transferred to 
part 1 to begin computation of the next step. 


The problem is terminated when the time limit TMAX is reached. This check 
is done in subroutine STEP. Had the problem exceeded the step number limit 
STEPMX, it would have terminated at that point. In either case, control is re- 
turned to MAIN 1 in segment 1 to print out the computing time and begin the next 
problem. When no data for another problem are given, the execution is terminated 
(i.e., control is returned to the monitor by subroutine INPUT as a result of an 
end of file on tape 7). The output of the last step is: 


STEP** 821. + 45. EC CENTRIC ITY* 
TIME= 42605.000 SEMILATUS R.« 
JDAY* 2440000.4927 MEAN ANOMALY* 
ALFA* 0. PATH ANGLE* 


2. 37578762E-04 OMEGA* 
6898571.50 TRU A* 
1.57042252 NODE* 
1 . 36122511E-02 INCL* 


1.57668670 

1.57089765 

0 . 

0 . 


V* 7599.09540 
VX* 43.7259269 
VY=s-7598 . 96967 
VZ=-0. 


R= 6898571.62 
X=-6898447 . 56 
Y*-41333 .9687 
Z*-0. 


REFER*EARTH ORBIT 1 
RMASS* 3846.70401 
REVS.* 7.50095356 
DELT* 263.055664 
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The time histories of several trajectory parameters for this example are 
shown as solid lines in figure 3. The oscillations of the eccentricity and mean 





Figure 3. - Time histories of several trajectory parameters for example 1. 


anomaly cause a rather small step size; as is noted in the figure. To indicate 
how exercising care in selecting the input can increase the computational effi- 
ciency; the same problem may again be run with the following initial values (ac- 
cording to ref. 9) of eccentricity and mean anomaly: 

2(F/m 0 )p 2 ^ n , e Q V 0 

e 0 - u ’> “0 - 2 " 3e 0 " 2Ig c 

The input cards for this case make use of the algebraic properties of the input 
routine to compute the desired value of these parameters. The cards are: 


$0ATA = 1 *$TABLE*41=RMASS*47=P, 5 =S IMP. 33* $$ IDENTIFICATION AND 

FLOW,10=DELT*30=TMAX*20=STEPMX«21=STEPS/ $$ TABLE DEFINITION 

RMASS=3850,$IMP=2540*rLOW=7.7361935E-5 $$ VEHICLE MASS* ISP* MASS FLOW 
P=6.86E6*7MAX=42605*STEPMX=200C $$ SEMILATUS-RECTUM* TIME LIMIT* STEP LIMIT 

DELT=1500*STEPS=200 $$ INITIAL STEP SIZE* OUTPUT EVERY 200TH STEP 

STABLE ,42=E ,46=MA/ E=2* 5 . 005 1 94 8E-4*P*P/ 3 . 98 3667E14 $$ ECCENTRICITY 

MA = -7620 *429/5 1 MP/9* 80665 -6 *E + 3« 141 5926/2 *STEPS = 5 55 MEAN ANOMALY *OUTPUT CONTROL 
SDATA=1* $$ LAST CARD 


The dashed lines in figure 3 show the time histories of the same trajectory 
parameters when initial values of e and M given immediately preceding are 
used. The increase in average step size is 20 to 1. To compare the accuracy of 
this approximation with the exact case (eQ = Mq =0); the final time was chosen 
when the corresponding orbit positions were identical (when the true anomalies 
were equal). At t = 42;605 seconds; the orbit positions are nearly identical; 
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and ; at this time, the values of position and velocity may he compared as fol- 
lows : 



Case A: 

Case B: 


e 0 = ^ = 0 

ep and Mq ^ 0 

Radius, m 

6898571.62 

6898571.56 

Velocity, m/sec 

7599.09540 

7599.09546 

Number of steps 

821 

39 


For most purposes the two answers would he accepted as equivalent and case B 
would he preferred because of the smaller computer time required. 

Example II; Lunar impact probe . This example of a lunar impact probe il- 
lustrates the use of the ephemeris tape and the control parameters needed to con- 
sider the effects of perturbing bodies, atmospheric forces, oblateness, rotating 
Earth, and thrust. No effort was made to optimize this trajectory but rather to 
use at least plausible values for illustrative purposes. 

Suppose the probe is launched at Cape Canaveral on December 7, 1961 by a 
three-stage vehicle with stage parameters as shown in the following table: 


Parameters 

Stage 

i 

2 

3 

4 

Initial mass, m q, kg 

150,000 

52,500 

23,625 

945 

2 

Engine exit area, A e , m 

3.0 

1.0 

.5 

(Coasting 

Vacuum specific impulse, I, sec 

300 

420 

420 

pay- 

.oad) 

Propellant loading, Wp/W q 

.65 

.55 

.96 



Propellant fraction, Wp^/Wp 

.9 

.9 

.91765873 



Propellant .flow rate, -ft, kg/sec 

750 

125 

56.25 



Burning time, t^ = Wp^/W, sec 

117 

207.9 

370 



Aerodynamic reference area, S, m^ 

7.5 

4.0 

2.0 

2 

.0 


Figure 4 shows the assumed variation of Cp q, Cp and Cp with Mach number as 
well as the angle-of -attack schedule. 

The vehicle will be flown as follows: First, a short nondrag vertical 

flight, after which the desired velocity orientation will be set, and then a turn 
determined by gravity and the angle-of -attack schedule until first-stage burnout. 
The second and third stages follow the same turn pattern. The final stage con- 
sists of the payload. The staging will be accomplished by treating each stage as 
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Mach number, Njyj 

Figure 4. * Angle-of -attack schedule and variation of drag and 
lift coefficients with Mach number. 


a single flight, with the turnout conditions of the previous stage used as ini- 
tial conditions. The chosen integration mode will he rectangular for the powered 
flight hut the mode of orhit elements will he used for the coast portion. Other 
bodies considered besides the Earth and the vehicle are the sun, the moon, and 
Jupiter. Jupiter is included to illustrate the use of ellipse ephemerides. The 
sun axid moon will illustrate the use of the tape ephemeris. 

The correct firing direction and launch time remain to he determined. This 
determination can he made by finding approximate values and then adjusting these 
values after one or more shots are fired. The adjustments could he made by an 
iteration scheme programmed internally to make a closed system. For this exam- 
ple, however, they were made by hand by firing several shots at various azimuth 
angles close to an estimate obtained by using reference 10 and an ephemeris. 

From a plot of the z-direction cosine of the vehicle-moon distance against 
vehicle-Earth distance, the azimuth angle that will intersect the moon orbit can 
be determined. The correct launch time is found by using the previously deter- 
mined azimuth angle and various times of day to determine the time of day at 
which the vehicle intersects the correct position in the moon orbit (location of 
the moon). This type of analysis gives an azimuth angle of about 78.9 and a 
time of day of about 7.94 h E.T. (E.T. is ephemeris time which is approximately 
equal to Greenwich mean time.) For the present purpose, these values will be 
used. 
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The problem begins by constructing the merged ephemeris tape for the sun and 
moon. This is done by subroutine TAPE in conjunction with the input shown as 
follows : 


$DATA=300»$TABLE»2=TAPE3,17=ELIST»29=TBEGIN»30=TEND/ $$ ID. AND TABLE DEFINITION 

TAPE3=0 $$ NECESSARY TO MAKE TAPE 

TBEGIN=2437640.5 $$ JULIAN BEGINNING DATE 

TEND=TBEG I N+5 $$ JULIAN ENDING DATE 

ELIST=(A3 (SUN, (A4JM00N $5 LIST OF DESIRED EPHEMERIS BODIES 


After the merged ephemeris tape is constructed, the clock is read, the standard 
data are initialized, and the first -stage input is loaded as shown: 


$DAta=1*$TAB,104 = LAT» 105 = LONG »106=AZI . 107=ELEV, 108* ALT* 
Z8.=IMODE,31=DTOFFJ,32=TOFFT»811=BODYCD»26=ATMN.29=RATM»459= 
ROTATE,41=RMASS,5=SIMP,33=FLOW.35=AREA,24=AEXIT.27=OBLATN»941= 
El IPS, 601=COEFN,238.=ICC, 37=ERE! r ,17=ERLIMT,19=CLEAR,30=TMAX,20= 
STEPMX.7=TKICK»10=DELT» 103.=MODOUT,23=DELMAX.22*TMIN»21=STEPS./ 


$3 STAGE 1 
$5 ID. AND 
$$ TABLE 
$$ DEFINITION 
$$ 


$$ 

LAT=28.280,LONG=-80.571 ,ELEV=89.7 $$ LATITUDE. LONGITUDE# ELEVATION 

AZI=78.9,ALT=10,IMODE=4 $$ AZIMUTH, ALTITUDE. INTEGRATION MODE 

DTOFFJ=24 37640. 5 »T0FFT=7. 94/24 $$ TAKE-OFF DATE AND FRACTION OF DAY 

B0DYCD=( A5IEARTH. (A4 ) MOON . (A6 1 JUPITE. ( A3 ) SUN *5 BODY NAMES. 1ST IS ORIGIN 
AT MN=(A5)EARTH.RATM=IEll.R0TATE=7. 29211585 E-5 $$ ATMOSPHERE NAME .RADIUS .ROTAT ION 
RMASS=150000 .SIMP =300. FLOW=750 $$ VEHICLE MASS. ISP (VAC) .MASS FLOW RATE 
AREA=7.5»AEXIT=3.0»0BLATN=(A5)EARTH $$ DRAG AREA, ENGINE EXIT AREA, OBLATE BODY 
EL IPS=(ALF6 ) JUPITE. ( ALF3 ) SUN, .9547861E-3 ,4. 81E+10, 5. 1913995. $$ ELLIPTIC DATA 

.0486288 , . 176 5935, .056971884, .405 87194 .24 33964. , . 6664 » 4333. 7 1 53$$ FOR JUPITER 
COEFN = 0, .4 ,0, .6, 1 , 1 . 15306 ,-.16326, .010204,8 » .5. ,» 100. ,10, » » $$ AERO. COEFF. AND 

100, , .025 , ,,100, ,,,,15, -.6, .04, ,40, .7, » , 1 17 , , , , 1E6 , ICC = 24 ♦ 14 , 19, 1 $$ INDICES 
EREF=lE-5 »ERL IMT = 5E-5 *CLEAR = 1 $$ REFERENCE ERROR, LIMIT ERROR .STDATA BY-PASS SWT 
TMAX=117,STEPMX=250 $$ MAXIMUM ALLOWED PROBLEM TIME AND STEP NUMBER 
TK1CK=10»DELT=2 $$ TIME OF THE VERTICAL NON-DRAG STEP. 1ST INTEGRATION STEP SIZE 
MODOUT=2 ,DELMAX=60, $$ MODE OF OUTPUT, TIME INTERVALS OF OUTPUT 


The value of MODE is set equal to 4, which causes execution of subroutine 
TUDES. TUDES transforms the spherical Earth coordinates into rectangular coordi- 
nates, which are the variables of integration. In addition, TUDES computes the 
closed-form solution for the initial vertical nondrag step. From this point on, 
the trajectory is integrated with the initial orientation specified by the spher- 
ical coordinates. The small error introduced by this procedure is offset by 
avoiding the complications associated with integrating the takeoff. One such 
difficulty is the thrust-direction specification when the velocity is zero, espe- 
cially if the origin body is rotating. 

Subroutine ORDER reorders the list of bodies putting the sun before Jupiter 
(i.e.j the sun's position relative to the vehicle must be found before Jupiter's 
relative position can be computed). The elliptic data for finding Jupiter's 
position are modified somewhat and relocated according to the computed body list. 
After calculating the gravitational constants, control is returned to MAIN 1. 
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The atmosphere belongs to the body at the origin (Earth) so that the rota- 
tion rate and atmospheric radius are set. The final duty of MAIN 1 is to posi- 
tion the merged ephemerides tape at the beginning of the correct ephemeris. In 
this case, only one merged ephemeris was constructed; nevertheless, it still must 
be identified and spaced to the beginning of the data. 

Control then passes to MAIN 2, where integration takes place in the same 
manner described in example I. Additional subroutines called from EQUATE are 
EPHMRS, ELIPSE , ICAO, AERO, THRUST, and OBLATE. Subroutine EPHMRS is responsible 
for computing the perturbations that result from bodies other than the origin 
body. This computation is accomplished by determining the perturbating body 
position through use of the merged ephemeris tape or subroutine ELLIPSE. 

The AERO subroutine determines the aerodynamic accelerations through use of 
quadratic equations for the lift and drag coefficients and subroutine ICAO, which 
determines density, pressure, and temperature as functions of altitude. Oblate- 
ness accelerations are found in subroutine OBLATE. The thrust direction is de- 
termined by subroutine THRUST, while the thrust magnitude is computed in EQUATE 
as ihg c I - PAg. 


The first vehicle stage integration is terminated by subroutine STEP when 
t = 117 seconds. Control is then transferred to MAIN 1, where the following in- 
put initiates the second vehicle stage integration: 

$D= 1 »RMASS = 52 500 »SIMP=420 »FLOW= 1 25 »TMAX=TMAX+207»9» AREA =4 »AEXIT=1 $$ STAGE 2 


Integration of the second stage proceeds in a manner similar to the integra- 
tion of the first stage and is terminated when t = 324.9 seconds. The third- 
stage data are similar to the second-stage data and are as follows: 

$D=1.RMASS=23625*FL0W=56.25»DELMAX=100*TMAX=TMAX+370»AREA=2.AEXIT=.5 $S STAGE 3 


The fourth stage differs from the preceding stages since the thrust is turned off 
and integration proceeds in orbit elements rather than in Cartesian coordinates. 
Output occurs every 6 hours until t = 1 day; then it occurs at every tenth step. 
Also, the error-control data are printed (therefore, make EREF negative). The 
fourth-stage input is as follows: 


$D=1*RMASS=945.0»DELT=3600*FLOW=0*TM AX =172800 S$ STAGE 4 

IMODE=— 2 » EREF=- ( ! $$ INTEGRATE ORBIT ELEMENTS* RECORD ERROR DATA. 

MODOUT=3*DELMAX=DELT*6.STEPS=10»TMIN=86400 $$ OUTPUT EVERY 6 HOURS UNTIL TIME = 

S$ 86400 * THEN EVERY 10TH STEP 


$D=1» S$ LAST CARD 


About l/2 day later the vehicle is close enough to the moon that the coordi- 
nate system origin is translated to the moon. This translation is accompanied by 
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a shift in integration mode to Cartesian coordinates, since the vehicle is ap- 
proaching the moon far out on a hyperbolic leg. The last step output is repro- 
duced as follows : 

STEP® 184. + 17. ECCENTRICITY® 10.6771772 0MEGA®-3 . 22087839 

TIME® 172800.00 SEMILATUS R.» 3.16835663E 09 TRU A® 1.16945998 

JDAY® 2437642.8306 MEAN AN0MALY®-18 . 8108633 NODE® 0.77242933 

ALFA® 0. PATH ANGLE® 62.2506247 INCL® 0.51408862 

MOON R® 2 . 556007 9E 08 -0.391661 0.734785 0.553798 
JUPITE R® 8 . 4571112E 11 0.581702 -0.741635 -0.334068 

V® 3938.07312 R= 6.12713230E 08 REFER=EARTH RECTAN 3 

VX= 2403.45856 X® 1.27258404E 08 RMASS® 944.999992 

VY=-2445. 76614 Y®-5 . 36514068E 08 REVS.® 0.78706574 

VZ=-1936. 50073 Z=-2. 67161870E 08 DELT® 5887.73633 

SUN R= 1 , 4668335E 11 -0.229169 -0.893118 0.387068 

At this time the vehicle is again primarily under the Earth's influence after 
missing the point mass moon by 1.2xl0 6 meters. 

Lewis Research Center 

National Aeronautics and Space Administration 
Cleveland, Ohio, September 6, 1962 
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APPENDIX A 


SYMBOLS 

— ► ^ — ► 

A relative angular momentum per unit mass, r x V r (appendix B) 

Ae engine exit area, m 2 

a i,j coefficients for quadratic functions 

Cp total drag coefficient 

Cp^Q zero angle- of -attack drag coefficient 

C'Dji induced drag coefficient 

C L lift coefficient 

D drag force, nevtons 

E eccentric anomaly, radians 

e eccentricity 

F thrust force, nevtons 

fl,fg functions of Mach number 

g c gravitational conversion factor, 9.80665 m/sec 2 (sometimes referred to 

as standard Earth gravity) 

h altitude above Earth r s surface, m 

I vacuum specific impulse, sec 

i orbit inclination to mean equator of 1950.0, radians 

J second harmonic coefficient in oblateness equations 

K fourth harmonic coefficient in oblateness equations 

k 2 universal gravitational constant, 1.32452139X10 2 ®, 

m 3 /(sec^) ( sun mass units) 

L lift force, ne'wtons 

M mean anomaly, radians 

m object mass, kg 
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m i 

m r 

% 

p 

, p 
p 

q 

R r 

r 


r i 

S 

T 

t 

U 

U x ,U y ,U z 

u 

V 

V' 


v 


W 


Wpf 

X 


x,y,z 


mss of 1 th perturbating body, sun mss units 
mss of reference body plus m, sun mss units 
Mach number 

atmospheric pressure, newtons/m 2 
f 1 xl (appendix B) 
power, w 

semilatus rectum, m 

dynamic pressure, ~ p(V') 2 ^ newtons/m 2 

radius of reference body, m 

radius from origin to object, m 

radius from origin to i*"* 1 perturbating body, m 

aerodynamic reference area, m 2 

t empe rat ur e , °K 

time, sec 

gravitational potential 

x,y,z accelerations . due to gravity, m/sec 2 

CD + V 

absolute velocity, m/sec 
relative velocity, m/sec 
true anomly, radians 
object weight, newtons 

propellant loading, fraction of mss that departs during a stage 

propellant fraction, fraction of W p used for propellant 

forces acting on object other than gravity, thrust, lift, drag, and 
perturbations due to perturbating bodies 

components of r, m 
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a 


angle between thrust and velocity vectors (sketch (a)), deg 


0 angle of rotation of thrust out of orbit plane (sketch(a)), 

T) power efficiency factor 

[x k^ m r 

p atmospheric density, kg/m 3 

cd argument of pericenter, radians 

cd origin body rotation rate, radians/sec 

& equatorial longitude of ascending node, radians 

Subscripts : 

0 initial value 

1,2, 3, 4 values at consecutive points along trajectory 



APPENDIX B 


VECTOR RESOLUTION 
Relative Velocity 

The relative velocity is defined as the velocity of the object with respect 
to the origin body. If the origin body is assumed to rotate about the z-axis, 
this velocity is given by 

V' = V - CD X r (Bl) 


In x,y,z component form, 

v x “ V x + ay ( B2a) 

V ' = V - ok ( B2b) 

«y «y 

V* = v z ( B2c) 

In the following sections, the atmosphere of the origin body is assumed to ro- 
tate as a solid body at the rate 


Thrust Resolution Along x,y,z Axes 

The thrust direction is specified with respect to the relative velocity 
vector V' by the angles a and J3, as shown in sketch (a). For resolution 
of thrust vector into x,y,z components^ it is convenient to define vectors A 
and P normal to and within the r, V' plane, respectively, such that V', 
A, and P form an orthogonal set. Thus, 


A = r X V' = 

Relative angular momentum per unit mass 

(B3) 


tc 

X 

ti> 

in 

tPn 

(B4) 

*-> — ► — ► 

The thrust vector can then "be resolved in the V 1 , A, P set as i 



F • V' = FV' cos a 

( B5a) 


F • A >= FA sin a sin £ 

( B5b) 


F • P = FP sin a cos p 

( B5c) 

Solving for F yields 



F = (V' cos a A X 

? + A sin a sin p ? x V + ? sin a cos (3 ?) 

(B6) 


P- 
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or, in x,y, z component form, 



cos a(AyP z - A z P y ) + A sin a sin 0(PyV^ - P z Vp 

+ P sin a cos 0 P x ] 

cos a(A z P x - A X P Z ) + A sin a sin 0(P z Vi - P X V Z ) 

+ P sin a cos 0 P y J 

cos a(A x Py - AyP x ) + A sin a sin P(P x Vy - PyVx) 


+ P sin a cos 



(B7a) 


(B7b) 


( B7c) 


Aerodynamic Lift and Drag Resolution Along x,y,z Axes 

— > 

The drag vector is alined with the relative velocity vector V' and is 
therefore given in x,y,z components as 


V 


V' 


D = -D £ - D £ - D ^ 


(B8) 


The lift vector L may be resolved into components along the previously 
defined orthogonal set V', A, and P by the following relations* 


L • V' =0 (B9a) 

L • A = LA sin 0 (B9b) 

L • P = LP cos 0 (B9c) 

Solving for L yields 

L = (A sin 0 P X V' + P cos 0 ?) (BIO) 

P 2 

or, in x,y,z component form, 

L x = ~[A sin 0(PyV^ - P z Vp + P cos 0 P x ] (Blla) 

Ly = ^[A sin 0(P z Vi - P X V^) + P cos 0 P y ] (Bllb) 

L z = ^[ A sin P( p x v y " P y V ^ + p cos 0 P z ] (Bile) 
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APPENDIX C 


TRANSFORMATION EQUATIONS BETWEEN RECTANGULAR 
COORDINATES AND ORBIT ELEMENTS 


Z 



From spherical trigonometry used in reference to the celestial sphere 
shovn in sketch (c) the following relations may he derived for the position 
coordinates : 


x = r(cos SI cos u - sin SI sin u cos i) 
y = r(sin SI cos u + cos SI sin u cos i) 
z = r(sin u sin i) 

where 


1 + e cos v 
u = O) + v 


and v is found from the relations 


and 


cos v 


cos E - e 
1 - e cos E 


M = E - e sin E 


(Cla) 

(cib) 

( Clc) 

(C2a) 

(C2b) 

(C2c) 

(C2d) 
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The velocity components may he obtained by differentiating the position equations 


using the two-body relations u = v 


- - Vkp 

n = v = » 


and 


* P 


e sin v: 



(C3a) 

(C3b) 

(C3c) 


where 

N = e cos co + cos u 
Q = e sin co + sin u 


( C4a) 
(C4b) 
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APPENDIX D 


RUNGE-KUTTA AND LOW-ORDER INTEGRATION SCHEMES WITH ERROR CONTROL 


The Runge-Kutta formula used is of fourth-order accuracy in step size h. 
It Is of the form 



x l = | (k, + 


2k c 


2k 3 + k 4 ) 


(Dl) 


where 


X = a dependent variable 


= increment in the dependent variable 


J 1 


hg = increment in the independent variable t 


ki = h^Ct^xp 


k 


2 


k 3 


h 2 X 2 1*1 + ~T> X l + 


h 2 X 2 (W + ^ X X + 

h 2 x 2 (n + h 2 > X 1 + 



A lower-order formula may be found by utilizing the three derivatives at 
t = *^l- y tg . If h-^ = t*j^ - 1 q and hg = tg - t-^ the following Lagran- 

gian interpolation formula gives the derivative at any time t Q < t < tg: 


^ ^ ” "^l^(t " ^g) •j. ^ “ ^2^ , ,v ^ “ ^p) 

A — Af ■ * - — ■ _L v 


hi(hi + h 2 ) A 1 h x h 2 + X 2 h 2 (h! + h 2 ) 

Integration of this equation from t-j_ to t 2 yields 


(D2) 


X 


I 


2 

1 



+ Shpx-L 




(D3) 
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The difference in the increments over the interval hg between the Runge-Kutta 
scheme and the_low-order scheme may be divided by a nominal value of the depend- 
ent variable X to obtain the relative error 8 2 - Thus, 


& 


2 



(D4) 


The error is expected to vary as approximately the fifth power of h, which 
leads to 


8 = Ah 5 


(D5a) 


(where A is a suitable coefficient) or in the logarithmic form 

log 8 = A' +5 log h (D5b) 

where 

A' = log A (D6a) 


Let it be assumed that A' will vary linearly with t, the variable of integra- 
tion. Then A' at a time corresponding to t3 can be found from A' at two 
previous points tp and tg as 

Aj = AJ + ^7 (t 3 - t 2 ) (DGb) 

and if hg = (tg - tg) and hg = (tg - tp 


A 3 “ A 2 + ( A 2 " A l^ 


^3 

h 2 


(D6c) 


and on this basis 8g would be predicted to be 

log &g = + 5 log h 3 


(D7) 


It is desired that 63 should approximate 5, the reference error; therefore. 


log h 3 = jr (log & - A£) 


(D8) 


Each dependent variable has an associated relative error and vould lead to com- 
putation of a different step size for each variable } hovever, the maximum rela- 
tive error of all variables may be selected for S. Obviously, inaccurate pre- 
dictions of step size can occur vhen the maximum relative error shifts from one 
variable to another or vhen any sudden change occurs . When a step size produces 
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an excessively large error (5 > S um-t -h) , a reduced step size must be used. It 
may be obtained from the reference error 5 as 


h 3 



(log 5 - ApJ 


(D9) 


S tarting the integration . - The Runge-Kutta scheme is simple to start, since 
integration from X n to X n +q requires no knowledge of X less than X n . 

Since the error control coefficient A has no value at t = 0, however, a pre- 
diction of the second step size is difficult. To overcome this difficulty, two 
equal size first steps may be made before checking the error. The A for the 
first step may be arbitrarily set equal to the A for the second step so that 
h 3 may be predicted. The low-order Integration scheme equation in this case 
becomes, with h 2 = hq, 


X 



(DIO) 


Failures . - Should two consecutive predictions of the same step fail to 
produce an error 6 less than §limit.> a return to the starting procedure will 
be made with a third prediction on step size, which is no larger than one-half 
of the second estimate. The step-size control described here will operate stably 
with nearly constant error per step only for a well-behaved function. For most 
problems it will repeat a step occasionally to reduce a large error, and on sharp 
corners it will restart. This action is not regarded as objectionable. The ob- 
jective Is to attain a desired level of accuracy with a minimum total number of 
steps . 


32 



APPENDIX E 


GLOSSARY OF VARIABLES 


Variable 

COMMON 

location 

Definition 

A 

562 

Total angular momentum per unit mass, m 2 /sec 

A (3) 

559-561 

x,y,z components of angular momentum per unit mass, 



m^/sec 

A1 

236 

Error control parameter defined by eq. (D6a) at tq 

A2 

237 

Error control parameter defined by eq. (D6a) at tg 

ACOEF 1 

265 


Interpolation polynomial coefficients ^ for variable 

ACOEF 2 

266 


step size (coefficient of Xq, Xp, Xg in 




eq. (D3)) 

ACOEF 3 

267 

J 


AJC (3) 

233-235 

Runge-Kutta coefficients; set in STDATA 

ALPHA 

564 

Angle between velocity and thrust vectors, positive 



when thrust vector is outward (sketch (a)) 

ALT 

463 or 108 

Vehicle altitude above an elliptic Earth, m 

AMASS (30) 

881-910 

Permanent list of body masses (sun mass units) in 



order of PNAME list; set in STDATA; masses from 
ELIPS data "begin at AMASS (21) 

ANGLES (4) 

104-107 

Same as LAT, LONG, AZI, and ELEV, respectively 

AREA 

35 

Effective area used to compute lift and drag forces 



in AERO, m 2 

ASQRD 

563 

Square of total angular momentum, A^, m^/sec^ 

ASYMFT 

543 

See table II 

ATMN 

26 

See table II 

AW (4) 

261-264 

Runge-Kutta coefficients; set in STDATA 

AZI 

106 

Azimuth angle, measured east from north at local 



meridian, input in deg 

BETA 

565 

Angle between velocity-thrust plane and orbit plane 



(sketch (a)) 
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Variable 

COMMON 

location 

Definition 

BEX (14) 

801-813 

List of error data 

B4ASS (8) 

417-424 

Body masses selected from AMASS list in sequence 
corresponding to EKAME list 

BNAME (8) 

402-409 

Ordered list of BCD body names 

BODY CD (8) 

811-818 

Original unordered list of BCD body names read from 
cards 

BODY L (8) 

801-809 

Auxiliary ordered list of BCD body names 

CD 

797 

Total drag coefficient per unit area, sec^/m 

CDI 

795 

Induced drag coefficient per unit area, sec^/m 

CEX (800) 

801-1600 

Common extension; common used in segment 1 but not 
needed in segment 2 and therefore saved on drum 2 
during execution of segment 2 

CF (126) 

276-401 

Coefficients from ephemerides tape used to determine 
positions of perturbing bodies 

CINCL 

495 

cos i 

CIRCUM 

541 

Circumferential component of total perturbative 
acceleration 

CHAMP 

246 

Smallest critical radius within which object lies 

CL 

796 

Lift coefficient per unit area, sec^/m 

CLEAR 

19 

See table II 

CLOCK 

3 

Contains reading of clock (to compute time used for 
particular problem) 

COEFN ( 190 ) 

601-790 

Storage array for coefficients used to compute 
ALPHA, CL, CDI, CD, or other parameters 

COMPA (3) 

537-539 

Components of total perturbative acceleration in 
x,y,z coordinate system 

CON (9) 

576-581 

Constants in the oblateness equations ; set in STDATA 

CONSTU 

18 

See table II 
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Variable 

COMMOE 

location 

Definition 

CONSU 

36 

See table II 

COS ALF 

575 

cos a 

COS BET 

599 

COS P 

COSTRU 

493 

COS Y 

COSV 

497 

COS U 

DE 

162 

e 

DEL 

255 

Used to control output in STEP 

DELMAX 

23 

See table II 

DELT 

10 

Step size, sec 

DIKCL 

165 

i 

EM 

161 

m 

DMA 

166 

M 

DECIDE 

164 

a 

DESITY 

460 

Atmospheric density, kg/m^ 

DQMEGA 

163 

d) 

DRAG (3) 

531-533 

x,y,z components of the drag acceleration 

DTOEF J 

31 

Julian date of takeoff 

E 

42 

e 

E2 

260 

Largest of relative errors between Runge-Kutta and 
Simpson rule integration methods defined by 
eq. (D4) 

EFMRS (7) 

410-416 

List BCD body names whose positions are to be deter- 
mined from ephemerides-tape data 

ELEV 

107 

Elevation angle, measured outward, deg 
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Variable 

COMMON 

location 

Definition 

ELIPS (120) 

941-1060 

Ellipse data for perturbing bodies, read from cards; 



for each body there are 15 pieces of data 



[NOTE: SUBROUTINE ORDER then converts 15 pieces of 



data into working set of 15] 

EPAR 

245 

>Je 2 - 1 

EREF 

37 

See table II 

EPTiTMT 

17 

See table II 

ERLOG 

259 

Natural logarithm of EREF 

ETOL 

25 

See table II 

EXMODE 

244 

Eccentricity (used when IMODE = 3) 

EMONE 

243 

e - 1 

FILE 

249 

See table II 

FLOW 

33 

Rate of propellant flow, kg/sec 

FORCE (3) 

525-527 

x,y,z components of acceleration due to thrust 

GASFAC 

458 

Defined in AERO; set in STRATA 

GEOH 

465 

Geopotential, m 

GK2M 

469 

Gravitational constant, ji, m^/sec^ 

GKM 

470 

Square root of GK2M 

H2 

256 

Value of DELT for previous step 

I BODY (8) 

425-432 

Defined in SUBROUTINE ORDER 

ICC (5) 

238-242 

See table II 

IMODE 

28 

See table II 

INCL 

45 

i, radians 

IED (3) 

791-793 

Index set in STRATA 
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Variable 

COMMON 

location 

Definition 

INDERR 

491 

Number of sets of error data; set in ERRORZ for use 
in MAIN 1 

KSUB 

254 

Index of Runge-Kutta subintervals 

LENGTH 

257 

See table II 

LAT 

104 

Geocentric latitude, positive northward, deg 

LONG 

105 

Longitude relative to Greenwich, positive eastward, 
deg 

ma 

46 

M 

MBODYS 

441 

Number of perturbing bodies 

MODOUT 

103 

See table II 

NBODYS 

489 

Total number of bodies,, excluding vehicle 

NDUMP (4) 

268-271 

See table II 

NEFMRS (8) 

433-440 

Defined in SUBROUTINE ORDER 

NODE 

44 

fl, radians 

NPONG (5) 

11-15 

See table II 

NSKEP (4) 

272-275 

See table II 

! 

j 

NSTART 

247 

Internal control in MAIN 2 and EQUATE 

OBLAT (3) 

534-536 

x,y,z components of oblateness acceleration 

OBLAT J 

38 

Oblateness coefficient of 2 nd - harmonic 

OBLAT K 

39 

Oblateness coefficient of 4^ harmonic 

OBLAT N 

27 

See table II 

OMEGA 

43 

gd, radians 

OLDDEL 

225 

Value of DELT for previous good step 

ORBELS (6) 

227-232 

Array of output variables, either rectangular or 
orbit elements 
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Variable 

COMMON 

location 

Definition 

P 

47 

p, m 

P (3) 

571-573 

Defined in eq. (B4) 

PAR (3) 

798-800 

Defined by equations in SUBROUTINE THRUST 

PMAGN 

574 

Defined in equation form by SUBROUTINE THRUST 

PRESS 

466 

Atmospheric pressure , mb 

PUSH 

34 

Thrust force , newtons 

PNAA 

— 

ALF list of body names 

PNAME (30) 

821-850 

Permanent list of body names made from PNAA list in 
SUBROUTINE ORDER; ELIPS names begin at PNAME(2l) 

PSI 

462 

Path angle , angle between path and local horizontal , 
deg 

QVAL 

794 

Defined in SUBROUTINE AERO 

QX (3) 

522-524 

x,y,z perturbative acceleration components due to 
perturbing bodies, m/sec^ 

R 

442 

Origin to object radius, m 

RADIAL 

540 

Radial component of total perturbative acceleration, 
positive outward, m/sec^ 

RATIO 

600 

Ratio of adjacent step sizes 

RATM 

29 

! Radius of atmosphere, m 

RATMOS 

248 

Set equal to RATM when ATMN equals reference body 
name (ENAME (l)) 

RB (3,8) 

200-223 

x,y, z components of distance from all bodies to 
object, m 

RBCRIT (8) 

450-457 

List of sphere-of-influence radii of all bodies in 
BNAME list, xn 

RCRIT (30) 

911-940 

Permanent list of sphere-of-influence radii corre- 
sponding to PNAME list of body names, mj radii 
from ELIPS data begin at RCRIT(2l) 

RECALL 

9 

| See table II 
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Variable 

COMMON 

location 

Definition 

REFER (30) 

851-880 

List of reference bodies corresponding to PNAME list; 
reference bodies from ELIPS data begin at REFER(2l) 

REVOLV 

250 

Rotation rate (rad/sec) of reference body when 
ATMN = BNAME (l) 

RESQRD 

40 

Square of Earth's equatorial radius, m^; used in 
SUBROUTINE OBLATE; set in STDATA 

REVS 

490 

Revolution counter , used only for output 

RMASS 

41 

m, kg 

ROTATE 

459 

Rotation rate of a reference body, radians/sec 

REEL (8) 

442-449 

Distances between bodies and object in order of 
BNAME list, m 

RSQRD 

567 

Radius squared of object to origin, m2 

SAVE 

8 

See table II 

SIMP 

5 

Specific impulse, I, sec 

SINALF 

569 

sin a 

SINHET 

568 

sin P 

SINCL 

494 

sin i 

SINTRU 

492 

sin v 

SINV 

496 

sin u 

SPD 

253 

Seconds per day; set in STRATA 

SQRDK 

468 

Gravitational constant 

m^/( sec^) ( sun mass units); set in STRATA; value of 

1.495X10 44 in/AU (equivalent to solar parallax of 
8.80008445 sec of arc) was used to convert units 
from 2.959122083X10" 4 

(AU)3/(mean solar day)2(sun mass units) to 
1.32452139x1020 m^/( sec2) ( sun mass units) 

STEPGO 

101 

See table II 

STEPNO 

102 

See table II 
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Variable 

COMMON 

location 

Definition 

STEPMX 

20 

See table II 

STEPS 

21 

See table II 

TAB (189) 

1301-1489 

Table array of input variables and their common 
storage assignment; used by SUBROUTINE INPUT; 
room for 94 variables 

TABLT 

252 

Time measured relative to DTOEFT, days 

TAPE 3 

2 

See table II 

TDATA (126) 

276-401 

Same as CF 

TDEL (7) 

592-598 

One-half of time spacing between two particular 
adjacent entries of like body name on ephemerides 
tape; read from tape for each body 

TEST 

1 

See table II 

TFILE 

16 

See table II 

TIM (?) 

585-591 

Time for set of ephemeris data; read from ephem- 
erides tape; one for each body 

TIME 

48 

Time, t, independent variable, sec 

TM 

467 

Temperature, °K, times ratio of molecular to actual 
molecular weight 

TMAX 

30 

See table II 

TMIN 

22 

See table II 

TOPFT 

32 

Fractional part of takeoff day (Julian), days 

TRSFER 

224 

See table II 

TRU 

483 

v, radians 

TTEST 

251 

See table II 

TTOL 

226 

Time tolerance within which problem time minus TMAX 
must lie to end problem 

V 

475 

Velocity of object relative to origin V, m/sec 


40 



Variable 

COMMON 

location 

Definition 

VA3M (3) 

477-479 

x,y,z components of VQ 

VEFM (3,8) 

498-521 

x,y, z components of object velocity relative to all 
various bodies, m/sec 

VEL 

109 

Initial velocity at input 

VQ 

480 

Velocity of object relative to atmosphere, m/sec 

VQSQBD 

481 

(VQ)^, m^/sec^ 

VMACH 

471 

Mach number of object, Mjj 

VSQRD 

476 

V^, m^/sec^ 

vx 

42 

x- component of V; also In COMMON location C(472), 
m/ sec 

VY 

43 

y-component of V; also in COMMON location C(473), 
m/ sec 

VZ 

44 

z-component of V; also in COMMON location C(474), 
m/ sec 

X 

45 

x- component of R, m 

X (15) 

131-145 

Working set of integration variables 

XDOT (15) 

161-175 

Array of Integration derivatives 

XIFT (3) 

528-530 

x,y,z components of lift acceleration, m/sec^ 

XINC (15) 

146-160 

Increments of integration variables per step 

XP (3,8) 

176-199 

x,y,z components of perturbing body positions rela- 
tive to origin 

XPKCM (15,2) 

41-70 

Two 15-variable arrays j second is integrated and 
first contains values of integration variables for 
last good step; see table V 

XPRIMB (15,2) 

71-100 

Least significant half of double precision Integra- 
tion variables corresponding to XPRIM 

XWHOLE (15) 

544-558 

Temporary storage for integration variables 

Y 

46 

y-component of R 
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APPENDIX E 


LEWIS RESEARCH CENTER EPHEMERIS 
General Description 

The ephemeris data initially available on magnetic tape for use on the 
IBM 704 computer were from the Themis code prepared by the Livermore Labora- 
tory, evidently from U. S. Naval Observatory data. Later, an ephemeris was 
obtained from the Jet Propulsion Laboratory assembled as a joint project of 
the Jet Propulsion Laboratory and the Space Technology Laboratory. These data 
are given relative to the mean vernal equinox and equator of 1950.0 and are 
tabulated with ephemeris time as the argument. 

An ephemeris was desired for certain uses in connection with the IBM 704 
computer that would be shorter than the original ephemeris tapes mentioned and 
would be as accurate as possible consistent with the length. A short investi- 
gation of the various possibilities led to adoption of fitted equations. In 
particular, fifth-order polynomials were simultaneously fitted to the position 
and velocities of a body at three points. This procedure provides continuity 
of position and velocity from one fit to the next, because the exterior points 
are common to adjacent fits. Polynomials were selected rather than another 
type of function, because they are easy to evaluate. Three separate polynomi- 
als are used for the x, y, and z coordinates, respectively. 


Procedure Used to Fit Data 


The process of computing the fitting equations is as follows: 

(l) A group of 50 sets of the components of planetary position was read 
into the machine memory for a single planet together with differences as they 
existed on the original magnetic tape. The differences were verified by compu- 
tation (in double precision because some data required it); and any errors were 
investigated, corrected, and verified. Published ephemeris data were adequate 
to correct all errors found. 


(2) The components of velocity v x , v y , and v z were computed and stored 
in the memory for each of the 50 positions by means of a numerical differenti- 
ation formula using ninth differences, namely, 


X = (T x - T_ 1 ) 


AI_ 1 + AI +1 AIII_ 1 + AIII +1 AV_ X + AV +1 


12 


60 


AVII_ 1 + AVII +1 AIX_ 1 + AIX +1 


280 


1260 


(FI) 
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(See ref. 11 pp. 42 and 99 for notation,) Double-precision arithmetic was used 
for differences, but velocities were tabulated with single precision. 

(3) Coefficients C, D, E, and F in the fifth-order polynomial 

X = X Q + X 0 (T - T q ) + C(T - T q ) 2 + D(T - T Q ) 3 + E(T - T Q ) 4 + F(T - T Q f (P2) 

and its derivative 

X = X 0 + 2C(T - T q ) + 3D(T - T Q ) Z + 4E(T - T Q ) 3 + 5F(T - T Q ) 4 (F3) 

were found to fit a first point (which was far enough from the beginning point 
to have all differences computed) and two equally spaced points for each com- 
ponent of position and velocity. (The initial spacing is not important, as will 
he seen later. ) Spacing is defined as the number of original data points fit- 
ted by one equation. Single-precision arithmetic was used. 

(4) The coefficients C, D, E, and F in step (3) were then used in equa- 
tions (F2) and (F3) to calculate components of all positions and velocities 
given in the original data and lying within the interval fitted. These values 
were checked with the original data. Radius R and velocity V were com- 
puted at the times tabulated in the original data. If any component of the 
position differed from the original data by more than RXLO"^ or if any ve- 
locity differed from the original by more than VX10~^, the fit was considered 
unsatisfactory. 

(5) If the fit were considered unsatisfactory, this fact was recorded; and 
the spacing was reduced by two data points. Steps 2 to 4 were then repeated. 

If the fit were considered satisfactory, this fact was recorded; and the spacing 
was increased by two spaces. Steps 2 to 4 were repeated. The largest satisfac- 
tory fit was identified when a certain spacing was satisfactory and the next 
larger fit was not satisfactory. 

(6) The coefficients that corresponded to the largest satisfactory fit were 
recorded on tape in binary mode as follows: 
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Word 

number 

Data 

Mode 

Definitions and/or units 

1 

Planet name 

BCD 

Six characters (first six) 

2 

Julian date 

Floating point 

Date of midpoint of fit, 





Julian date 

3 

Delta T 



Number of days on each 





side of midpoint 

4 

F x 



a AU/dayE 

5 




a AU/day 4 

6 

-^X 



a AU/day 3 

7 

C x 



a AU/day 2 

8 

w 

X 



a AU/day 

9 

X 



a AU 

10 

F 

y 



a AU/dayE 

11 

Ey 



a AU/day 4 

12 

D y 



a AU/day 3 

13 

c y 



a AU/day 2 

14 

y 



a AU/day 

15 

y 



a AU 

16 

F 

z 



a AU/day 5 

17 

E z 



a AU/ day 4 

18 

D z 



a AU/day 3 

19 

c z 



a AU/day 2 

20 

z 



a AU/ day 

21 

z 

1 


a AU 


a Except for moon data, which are in Earth radii and days. 


(7) As soon as a set of coefficients was selected for an interval, addi- 
tional data were read from the source ephemeris tape and used to replace the 
points already fitted (except the last point). These data were processed as 
described in steps 1 and 2 so that the next 50 points were ready to he fitted. 
Steps 3 to 6 were then used to find the next set of coefficients, and steps 1 
to 6 were repeated until all data for all planets and so forth, were fitted. 


Data Treated 

The preceding process was applied to all data available at the time. For 
the moon, the technique usually led to the use of every point in the fitted 
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interval (i.e., only three points vere fitted). Thus, a check of accuracy was 
not available. The error on the attempt to fit the next greater interval 
(five points) was not excessive, however, and it is judged that the accuracy ob- 
tained from these fits is about equal to that held on the other bodies. 


Merged Ephemeris Tape 


Once all the positions and velocities of all the bodies then available were 
fitted, the coefficients were merged in order of the starting date of each fit. 
The resulting tape was written in binary mode with 12 sets of fits per record. 


The detail of this record is as follows: 


Set 1 


Set 2 


1 word: 
2 nci word: 
3^ word: 
4 th word: 


23 r< ^ word: 


FORTRAN compatible 

file number, fixed point in decrement 

planet name, code in BCD, first six characters 

Julian date, floating point 

etc. , according to list in paragraph 6 

21 words 

z 


24 th vor(i: planet name, code in BCD, first six characters 

25 word: Julian date, floating point 


44 th word: z 


Successive sets follow one another with a total of 12 sets. 


Set 12 
(last set) 


234 uu word: 
235 th word: 
254^ word: 
255^ h word: 
256 th word: 
End-of-record 


planet name 
Julian date, 
z 

zero 

zero 

gap 


floating point 


One record contains 256 words, the first is for FORTRAN compatibility, the 
second is a file number used for identification in the system. It Is a fixed 
point 2. The third is the beginning of the first set of data, and 12 sets fol- 
low, each with 21 words. The last word is the 256 th word (counting the FORTRAN 
compatible word) followed by an end-of-record gap. The remaining records are 
compiled in the same manner with an end-of-file recorded as a terminating mark. 

Because of the merging operation, all bodies are given In one list in a 
random order according to the starting date of the interval. The starting date 
is the Julian day (word 2) minus the half interval (word 3) (see procedure, 
paragraph 6). The entire ephemeris occupies about one-seventh reel of tape. A 
summary of data is given in table VII. 
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APPENDIX G 


INPUT -DATA REQUIREMENTS 

The procedure needed to run actual problems vith the aid of this routine is 
described herein. It is intended to permit a person vith a specific problem in 
mind to make a complete list of data required and to select desirable operating 
alternatives from those available to him. The details of this procedure are 
contained in the folloving instructions: 

(1) Provision has been made for two types of ephemeris data to specify the 
locations of celestial bodies that perturb the vehicle. They are ellipse data 
and ephemeris-tape data. If the problem does not Involve perturbing bodies (ex- 
cept a reference body) or if elliptic data are used for all the perturbing 
bodies, skip to instruction 5. 

(2) If the perturbing -body data are to be taken from an ephemeris tape, 
list the names of the ephemerides and Julian dates to be covered along vith the 
folloving auxiliary information: 

1 st card: $DATA = 300, STABLE, 2 = TAPE 3, 17 = ELIST, 29 = TBEGIN, 

30 = TEND/ 

Other cards: TAPE 3 = 0 

TBEGIN = ephemeris beginning Julian date 
TEND = ephemeris ending Julian date 

ELIST = (names of perturbing bodies in "ALF” format, see 
example in text) 

The ephemerides of all planets except Earth bear the name of the planet. The 
ephemeris giving the distance from Earth to the sun is called "sun," as Is 
astronomical practice. 

(3) If successive files on the ephemeris tape are to be made, punch the 
corresponding sets as follovs: 

$DATA = 300, TAPE 3=0, TBEGIN = , TEND = , ELIST = 

As many similar sets as are needed may be appended. 

(4) If ellipse data are to be loaded from cards, they are prepared later 
under instruction 12. 

(5) On the first execution after loading the routine, the common area Is 
cleared vhether an ephemeris tape is constructed or not. It is nov necessary 
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to load a table of variable names. Once loaded, this table will not be cleared 
again (except if the control variable TAPE 3 is set to zero). These names are 
for use on the input cards. If a different name is desirable for any variable, 
it may be changed in the table and where it appears on the input card (ref. 7). 
The cards are: 


SDAT A=1 »$TABLE*104=LAT* 10 5 = LONG »106=A2I*107=ElEV» 108-ALT*109=VEL»7=TKICK 
* 2 8 « = I MODE »45=X»46=Y*47=Z*42=VX»43=VY*44=VZ*42=E* 43=0MEGA * 44=N0DES * A 5 = 
INCLf46=MA,47=P*41=RMASS*31=DTOFFJ,32=TOFFT,48=TIME»811=BODYCD.16=TFILE» 
94l=ELlPS*27=OBLATN»38=OBLATJ » 39 = 0BLATK # 3 4 = PUSH » 5 = S I MP »33=F10W»24=AEXIT * 
56 5 = BETA t 601=C0EFN*2 38.= I CC » 26 = ATMN * 29 =R A TM »459=R0TATE * 35 = AREA * 37 = EREF * 
X7=ERLIMT*103.=MODOUT*30=TMAX»20=STEPMX»23=DELMAX»21=STEPS»22=TMIN» 1= 
TEST»268.=NDUMP*272.=NSKIP.257.=LENGTH,19=CLE/'R*8=SAVE,9=RECALL*10=DELT/ 


(6) The initial position and velocity of the vehicle may be given in any 
one of three coordinate systems. If the initial data are given in orbit ele- 
ments, skip to instruction 8. If the Initial data are given in rectangular co- 
ordinates, skip to instruction 7. If the initial data are given in Earth- 
centered spherical coordinates, the following variables should be punched: 

LA.T = latitude, deg, positive north of equator 

LONG = longitude, relative to Greenwich, deg 

ALT ss altitude above sea level, m 

AZI = azimuth angle, east from north, deg 

ELEV = elevation angle, horizontal to path, deg 

VEL = initial velocity, m/sec 

TKICK = size of initial vertical, nondrag step to facilitate starting, 
sec 

MODE = 4 

These geocentric coordinates are converted by subroutine TUDES to rectangular 
coordinates and MODE will be changed to 2 with its original sign. Skip to in- 
struction 9. 

(7) If the initial data are in rectangular coordinates, set the following 
variables : 

X = x-component of position in x, y, z coordinate system, m 
Y = y-component of position in x, y, z coordinate system, m 
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Z = z-component of position In x, y, z coordinate system, m 
VX = x-component of Telocity in x, y, z coordinate system, m/sec 

VY = y-component of velocity in x, y, z coordinate system, m/sec 

VZ = z-component of velocity in x,y, z coordinate system, m/sec 

IMODE = 2 

Skip to instruction 9* 

(8) If the initial data are in orbit elements, set the following variables: 

E = eccentricity 

OMDSA = argument of pericenter, radians 

NODES = longitude of ascending node (to mean vernal equinox of 1950.0), 
radians 

INCL = orbit inclination to mean equator of 1950.0, radians 
MA = mean anomaly, radians 
P = semilatus rectum, m 
IMODE = 1 

(9) Integration is performed on either rectangular variables or orbit ele- 
ments. If the Initial data are of the same type as the desired Integration 
variables, the positive sign on IMODE, as given in instruction 8, will signal a 
matching condition; but if the desired Integration variables are of the opposite 
type to the Input variables, a minus sign should be affixed to the value of 
IMODE. Note that in the case of geocentric coordinates, an automatic conversion 
to rectangular coordinates Is effected. To convert geocentric coordinates to 
orbit elements requires IMODE = - 4, whereupon subroutine TUDES will convert the 
geocentric coordinates to rectangular coordinates, IMODE will be set to -2, and 
then In MAIN 2 the further conversion to orbit elements will be sensed with 
IMODE finally being set to +1 by the program. 

(10) To specify vehicle mass and takeoff time, set the following variables: 
RMASS = mass of vehicle, kg 

DTOFEJ = Julian day number 
TOEFT = fraction of day 

TIME = time from previously set Julian date, sec 
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Takeoff occurs at the instant corresponding to the sum of the last three quan- 
tities. If a specific date or time is not required, these variables may be 
skipped. In that case, the STDATA subroutine sets DTOFFJ to 2440 000. 

(11) To specify the origin and any perturbing bodies, list them as BODYCD 
(list of body names in "AIF If format, see text example). The first body In the 
list Is taken to be the reference body. The distances between the bodies in 
this list must be computable from either ellipse data (instruction 12) or 
ephemeris -tape data (instruction 2). There may be no more than eight names In 
the list. Also, If the ephemeris tape is being used, the correct file must be 
found on It. For this purpose, set TFILE = desired ephemeris tape file. The 
ephemeris files were numbered in sequence when written in instruction 2. If 
TFILE is not given, it will be set equal to 1.0 by the STDATA subroutine. 

( 12 ) For each body whose path is represented by an ellipse, a 15-element 
set of data must be loaded. A 15-element set consists of: 


1. body name In ,f AIF ff format (maximum of six characters) 

2. reference body name In fr AIF" format (maximum of six characters) 

3. mass of body, sun mass units 

4. radius of sphere of influence, m 

5. semilatus rectum, AU 

6. eccentricity 

7. argument of pericenter, radians 

8. longitude of ascending node (to mean vernal equinox of 1950.0), 

radians 

9. orbit inclination (to mean equator of 1950.0), radians 

10. Julian day at perihelion 

11. fraction of day at perihelion 

12. period, mean solar days 



15. j 

It is convenient to punch a 15-element set in sequence and to separate the 
elements by commas on as many cards as are required. Several sets may then be 
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loaded consecutively. The order of the sets is immaterial. Ellipse data, if 
present, take precedence over ephemeris-tape data. The sets are loaded consecu- 
tively, in any order, as follows: 

ELIPS = set 1, set 2, set 3, . . ., set nj n < 8 (see example in text) 

(13) To specify the initial integration step size, set 
DELT = initial integration step size 

If no value of DELT is given. It will he set to TMAX/lOO by MAIN 1. 

(14) If oblateness effects of the Earth are to be included, set 
OBLATN = (ALF5) EARTH 

(15) If thrust forces are present, set either 

(1) PUSH = thrust magnitude, newtons (for m = 0) 
or 

(2) SIMP = specific impulse (vacuum), sec 
FLOW = mass-flow rate, m, kg/sec 

For either choice, set 

2 

AEXIT = engine exit area, m 

Also, the thrust orientation must be specified by setting 
BETA = angle (3, deg (see sketch (a)) 

COEFN (i) = angle-of-attack schedule, a = a(t) (see instruction 17) 

ICC = fixed-point integer (see instruction 17) 

For the special case of tangential thrust, none of the last three variables need 
be set. 

(16) If aerodynamic forces are present, set 

ATMN = name of body that has atmosphere, in TT AIF" format 

RATM = radius above which atmospheric forces are not to be consid- 
ered, m 

ROTATE = atmospheric-rotation rate, radians/sec (7.29211585X10 ^ for 
Earth) 
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AREA = reference area, m 2 

BETA = angle J3, deg (see sketch (a)) 

COEFN (i) = angle-of attack schedule, a = a(t), Cp/sin a, Cp q, and 
Cp 1 /Cl curves (see instruction 17) 


ICG = fixed-point integers (see instruction 17) 


(17) If neither thrust nor aerodynamic forces are present, skip to instruc- 
tion 18. The relations a,(t), Cp/sin a, Cp Q , and Cp -j ,/c|; are assumed to he 
quadratic functions that involve coefficients which are located in the COEFN(J) 
array. The arrangement of these coefficients is best explained by an example. 
Suppose the functions a(t) is as follows: 


r , ,2 

a ll + a 12 t + a 13 t 


a 21 + a 22 t + a 23 t 


a = ' a 31 + a 32 t + a 33 t2 


(ti < t < t 2 ) 
(t 2 < t < t 3 ) 


(t 3 < t < t 4 ) 


etc. 


etc. 


The coefficients a i,j should then be loaded into the COEFN(J) array as: 


COEFN(j) - tp, a u , a 12 , a 13 , t 2 , a 21 , a 22 , a 23 , t 3 , a 31 , a 32 , a 33 , t 4 , . . . , t n 


Furthermore, additional sets of coefficients for the other functions may simply 
he added to the COEFN(j) array, which results in a string of sets of coeffi- 
cients, and can he represented, for example, as; 

COEFN(j) = a coefficients, CjJ sin a coefficients, Cp q coefficients, 
etc* 9 

= tp, a 11? a 12 . . ., t n , N^p, bpp, bp 2 , . . ., %^ k , etc. 

The starting point in the COEFN(j) array of each function must also be loaded to 
identify the correct region of coefficients. To this end, the following array 
must also be loaded: 

ICC(l) = fixed-point value of J where a coefficients begin 
ICC(2) = fixed -point value of J where Cp/sin a coefficients begin 
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ICC (3) = fixed-point value of J where C^ coefficients begin 

ICC (4) = fixed -point value of J where C^ Q coefficients begin 

For this purpose, all values In the COEFN(j) array are called coefficients (i.e. > 
the t's and the %'s are coefficients). The sequence of the sets is arbi- 
trary, since changing the sequence requires only a change in the ICC(l) array. 

See Example II - Lunar impact probe section. 

(18) The size of the Integration steps Is determined primarily by the error 
control variables. These are loaded as: 

EKEF = error reference value; 6 in appendix D 

ERLIMT = maximum value of 6 that is acceptable on any particular step 

EREF Is always treated as a positive number; however, if it is loaded with a 
minus sign, this will cause error Information to be printed at the completion of 
the problem. If no error control data is loaded, STDATA subroutine will set 

EREF = lXlO -6 , ERLIMT = 3xl0“ 6 . 

(19) The output control offers a choice on the frequency of output data as 
follows : 

If MODOUT = 1, output will occur every n th step (n = STEPS) until 
t = TMIN, and then MODOUT is set equal to 2 by the program 

If MODOUT = 2, output occurs at equal time intervals of DEIMAX until 
t = TMAX 

If MODOUT = 3, output occurs at equal time intervals of DEIMAX until 
t = TWIN, then MODOUT is set equal to 4 by the program 

If MODOUT = 4, output occurs every step (n = STEPS) until 
t = TMAX 

TMAX — maximum time limit before problem is completed 
STEPMX = maximum step limit before problem is completed 
DEIMAX = time Interval between outputs 
STEPS = number of steps between outputs 
TMIN = time when MODOUT changes 

Note that output control may, at times, strongly influence the integration step 
size especially if MODOUT is 2 or 3 and DEIMAX is small. TMAX must be loaded. 

All others may be skipped; if so, STDATA will put MODOUT = 4, and STEPS = 1. 
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(20) For debugging operations and for occasional supplementary output, it 
may be desirable to obtain G-type format dumps. These may be obtained through 
strategic positioning of the FORTRAN calling statements CALL DUMP (ID, DATA, 
LENGTH) where ID is the identification number to appear in the output, DATA is 
the starting location of the dump area, and LENGTH is the number of consecutive 
words to be dumped. To actually obtain dumps at execution time, set 

TEST = total desired number of dumps 

NDUMP(j) = identification numbers of desired dumps, corresponding to 
ID ! s of calling statement, J < 4 

NSKIP(j) - number of skips to occur between dumps, NSKIP(j) acts upon 
NDUMP(J), J < 4 

LENGTH = number of consecutive words to be dumped 

Note NDUMP(j) will occur the NSKIP(j)^ time control passes through the calling 
statement and will occur every NSKIP(j) th time thereafter. If NSKIP(j) is omit- 
ted, it is taken to be 1. DATA may be a common location or the name of a rela- 
tive variable. If the value of a word to be dumped is zero, it is skipped. 

(21) For certain problems, it is desirable to save the initial data read in 
on cards or the data generated at the completion of a part of a problem. The 
saved data may then be recalled at a later time to be used as intial conditions 
for another problem. To prevent the "standard data" set from being loaded (and 
the accompanying common clearing loop), set 

$DATA =99, CLEAR = any nonzero number 

To save the initial data before the input is read in (i.e., the result of a pre- 
vious calculation), set 

SAVE = 2 

To save the initial data after the input is read in, set 
SAVE = 1 

To recall the saved data, set 

RECALL = any nonzero number 
CLEAR = any nonzero number 

By taking advantage of the place in the program where each discrimination is 
made, several useful combinations of these controls are possible (see fig. 2). 
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(22) When a transfer of origin occurs, provision has been made to read 
input into the program. This is done with the aid of $DATA - 101, followed hy 
the data statements desired. 

(23) Following is an input check list that may he helpful at execution 

time: 


INPUT CHECK LIST 8, 


Time and Mass 

Position and velocity 

(completely fill in one and only one block) 

Reference and perturbing bodies 


Rectangular 

Orbit elements 

Spherical 

BODYCD = 

Tape 

Elliptic 

DTOFFJ = 
TOFFT = 
CELT = 
TIME = 
RMASS = 

CO 

II 

W 

II tl II Q 
***£££§ 

E = 

OMEGA = 
NODES = 
INCL - 
MA = 

P = 

IMODE = 1 

i 

LAT = 

LONG = 
AZI = 

ELEV = 

AIT = 

VEL = 
IMODE = 4 

i 

TAPE 3 = 0 (b) 

TBEGIN = 

TEND = 

ELIST = 

TFILE = 

~ T 1 

ELIPS = 


Output control 


TMAX - 
TMIN - 
M0D0UT = 
STEPS = 
DELMAX = 
STEFMX = 


(c) 


Error control 


EREF = 
ERLIMT 


Restart 


SAVE « 
RECALL 
CLEAR = 


Thrust (d) 


SIMP = 
FLOW = 
PUSH = 


Atmosphere 


ATMN = 
RATM = 
ROTATE = 
AREA = 


COEF = 
ICC = 
BETA = 


Ohlateness 


OBLATN = 


Dump 


NDUMP = 
NSKIP = 
LENGTH = 
TEST = 


a The following standard data are loaded hy subroutine STDATA: 

nrOFFJ = 2440 000.0 MODOUT = 4 EREF = 1X10" 6 

IMODE = 1 STEPS = 1.0 ERLIMT = 3X10 “ b 

BODYCD(l) = (ALF5) EARTH STEPM5C = 100.0 TFILE = 1.0 

EMASS =1.0 

input 300, setting TAPE 3 = 0 is necessary to make an ephemeris tape. 
(^A value for TMAX is always required. 

( d )use either SIMP = and FLOW = or PUSH = 
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APPENDIX H 


PROGRAM LISTING 


C MAIN t — FLOW CONTROL PROGRAM FOR SEGMENT 1. THIS ROUTINE IS ENTERED FOR 

C EITHER ID THE START OF A PROBLEM OR 12) AN IN-FLIGHT ORIGIN BODY CHANGE. 

C THE REGION OF COMMON FROM 801 TO 1600, CEX, 15 NOT NEEDED IN SEGMENT 2 AND 

C IS THEREFORE SAVED ON DRUM 2 DURING EXECUTION OF THAT SEGMENT TO OPEN UP 

C THE ADDITIONAL 800 STORES. 

C 

COMMON C 


C 


C 


c 

c 

c 

c 


c 

c 

c 

c 


c 

c 


c 

c 


c 

c 

c 

c 


c 

c 

c 


c 

c 

c 

c 


c 


DIMENSION 

1 C 11600), BNAME 18), 

2 CEX 1800), BODYCD (8), 


NPONG 15), 
BEX (14) 


EQUIVALENCE 

1( CLEAR » C ( 19)), I CLOCK, C( 5)) 
21 TTOL , C I 226 ) ) , I IMODE,C( 28J) 
3 ( DEL,C(255) ) ,( £REF,C( 37)1 

4( TAB ,C I 1301 ) ) , IDELHAX ,C( 23)) 
5IREVOLV.CI250) 1 ,(ATM N ,C( 26)1 
6(R0TATE»C(459) ) , t NPONG, Cl ID) 
71 TF ILE » C { 16)), ( FILE ,C 1 249 ) ) 

81 cex,c(8om 


,( BEX, Cl B01) ) , I INDERR » C < 491) ) 
, (RECALL, CC 9)1, { SAVE, Cl 6)) 
, I LENGTH, Cl 257) ) , I TMAX,C< 30)) 
,( DELT.CI 10)), I ERLOGjCI 259)1 
, ( R ATM ,C( 29) ) , <RATMOS,CI240D 
, (BODYCD, C(811M, I BNAME, CI402)) 
,1 T APE3, C ( 2)), (TR5FER»C(224)J 


THE COMMON EXTENSION, CEX, IS RESTORED (JUNK IS BROUGHT IN UPON THE FIRST 
ENTRY). WHEN TAPE3-0.0, SUBROUTINE TAPE IS CALLED TO COMPILE THE 
EPHEMERIDES. SUBROUTINE TAPE ALWAYS SETS TAPE3-3. 

READ DRUM 2,0, CEX 
IF { TAPE 3) 2,1,2 

1 CALL TAPE 

WHEN AN IN-FLIGHT ORIGIN TRANSFER OCCURS, SEGMENT 1 IS CALLED WITH TRSFER 
•1.0. HERE, AN INPUT IS ALLOWED AND THEN CONTROL IS SENT TO REORDER THE 
BODY LIST. 

2 IF (TRSFER) 4,4,3 

3 TR5FER * 0. 

CALL INPUT! 101, C, TAB) 

GO TO 28 

PRINT OUT THE ERROR INFORMATION IF EREF HAS A - SIGN. 

4 IF (EREF) 5,10,10 

5 WRITE OUTPUT TAPE 6,8 
REWIND 4 

DO 6 1*1, INDERR 

READ TAPE 4, BEX 

6 WRITE OUTPUT TAPE 6 , 9 , I BEX ( J ) , J = l , 1 4 ) 

7 REWIND 4 
INDERR * 0 

READ DRUM 2,786,BEX 

8 FORMAT! 7H l STEP,6X,4HT1ME,6X, 4HDELT,7X, 2HA2 , 8X , 2HE2, 7X , 4HMASS, 6X, 

14HE ,VX,4X,8H0MEGA,VY,2X, 8HN0D6S *YZ» 3X,6HlNCL«X,5X, 4HMA, Y, 6X » 3HP , Z , 

24X , 1HK // ) 

9 FORMAT { F5. , IH + F3. , IP 1 1G10.2 , 12 ) 

PRINT OUT THE COMPUTATION TIME ELAPSED SINCE THE LAST ENTRY TO MAIN 1. 

10 CALL TIME1 (CLOCKl) 

IF (CLOCK) 11,13,11 

11 TUSED * CLOCKl - CLOCK 

WRITE OUTPUT TAPE 6 ,12, TUSED 

12 FORMAT! 15H0MINUTES USED -F7.1/1HU 

13 CLOCK - CLOCKL 

CALL IN THE STANDARD DATA IF CLEAR=0 . INPUT 99 IS BASICALLY AN AUXIALLARY 

input to allow a change in clear, if save-2. o, the data from common 

5 TO 115 IS SAVED. 

CALL INPUT (99, C, TAB) 

IF (CLEAR) 15,14,15 

14 CALI STOATA 

15 IF (SAVE-2.) 18,16,18 

16 DO 17 J-5,115 

17 CtJ+1485) » C(J) 

WHEN RECALL DOES NDT EQUAL ZERO, THE INITIAL DATA PREVIOUSLY STDRED BY A 
SAVE STATEMENT WILL BE RECALLED IN ORDER TO RESTART WITH THE SAME DATA. 

18 IF (RECALL > 19,21,19 

19 DO 20 J-5 , 115 

20 C1J) - CU+1485) 

INPUT l IS THE MAIN INPUT STATEMENT, DATA REAO IN HERE OVERWRITES ANY 
STANDARD VALUES SET BY STDATA. IF SAVE-1.0, THE INITIAL SET OF DATA FROM 
COMMON 5 TO US WILL BE SAVEO FOR LATER USE. 

21 CALL INPUT (1,C,TAB) 

IF ( SAVE-l. ) 24,22,24 

22 DO 23 J* 5 , 115 

23 C ( J+1485 ) - CIJI 

24 IF (DELT) 26,25,26 

25 DELT - T MAX/100. 

26 ERLOG - LOGF(ABSFIEREF) ) 

DEL - DELMAX 

TTOL * 5 E-8 *THAX 
BNAME ( l ) - BODYCD ( 1 ) 
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oooooonnononnrinnorioriooriooonnoonrinono 


C IF IMODfc IS 4, THE INITIAL DATA IS EARTH CENTERED SPHERICAL 

C COORDINATES ANO IS TO BE TRANSFORMED INTO RECTANGULAR COORDINATES. THIS 

C IS DONE BY SUBROUTINE TUOES WHERE* ALSO, AN INITIAL STEP MAY BE TAKEN TO 

C FACILITATE STARTING. THE COMMON EXTENSION, CEX, IS SAVED. SUBROUTINE 

C ORDER IS CALLED TO ORDER THE LIST OF BODIES, COMPUTE THE GRAVITATIONAL 

C CONSTANT, AND MODIFY ANY ELLIPTIC EPHEMERIS DATA. 

IF { XABSF I I MODE ) -4 ) 28,27,20 

27 CALL TUOES 

IMODE = XSIGNFI2, I MODE ) 

28 WRITE DRUM 2,0, CEX 
CALL ORDER 

CALL DUMP (I, C, LENGTH) 

C 

C IF ORIGIN BODY HAS AN ATMOSPHERE, SET ROTATION RATE AND ATMOSPHERE RADIUS. 

REVDLV » 0. 

RATMOS * 0. 

IF 1ATMN - BNAME(II) 31,29,31 

29 RE VOL V * ROTATE 
RATMOS * RATM 

C 

C POSITION THE EPHEMERI DE S TAPE AT THE BEGINNING OF THE CORRECT EPHEMERIS 

C BY MATCHING THE EPHEMERIS NUMBER READ FROM TAPE (FILE) WITH THE DESIRED 

C EPHEMERIS NUMBER (TFILE). THEN CALL IN SEGMENT 2. 

31 IF (FILE) 36,36,32 

32 CALL BKF I L E < 3 > 

33 READ TAPE 3, FILE 

IF (FILE-TFILEJ 34,36,32 
S 34 RTB 3 
S 35 CP Y 
S TRA *35 

S TRA *33 

S TRA *34 

36 CALL PONG ( NPONG ( l ) ) 

C 

C END OF THE FDRTRAN STATEMENTS. •»**#*** 


SUBROUTINE TAPE 

SUBROUTINE TAPE USES THE MASTER MERGED EPHEMERIDES TAPE (TAPE 8 AT LEWIS) 
TO COMPILE A WORKING EPHEMERIS TAPE l TAPE 3 AT LEWIS) WHICH CONTAINS ONLY 
THAT DATA NEEDED AT EXECUTION TIME. THIS MINIMIZES TAPE HANDLING DURING 
EXECUTION. THERE ARE 2 FILES ON TAPE 8, FIRST FILE HAS THE DATA AND IS 
IDENTIFIED BY THE SECOND WORD OF EACH 256 WORD RECORD (FIRST WORD IS THE 
DUMMY FORTRAN COMPATIBLE WORD, SECOND WORD-2). THE SECOND FILE IS ONLY 2 
WORDS LONG, FIRST WORD IS FORTRAN COMPATIBLE, SECOND WORD-3). 

MASTER FILE 1 — PLANETS (EXCEPT MERCURY ANO EARTH), SUN, MOON, AND 

EARTH-MOON BARYCENTER FROM SEPT. 25, 1960 TO ABOUT 2000. 
EACH EPHEMERIS COMPILED REQUIRES A SET OF INPUT 300 DATA. THE FIRST PIECE 
OF DATA WRITTEN ON A FILE IS THE FILE IDENTIFICATION NUMBER, FILE. EACH 
FILE IS NUMBERED CONSECUTIVELY STARTING WITH FILE-l. SINCE MOON DATA IS IN 
TERMS OF EARTH RADII, THE CONVERSION OF MOON DATA TO A.U. IS MADE BEFORE 
WRITING ON TAPE 3. THE COMMON USED IN SUBROUTINE TAPE IS LOCAL ANO ALL 
BUT T APE 3 IS CLEARED BY A FINAL CLEARING LOOP. 

FUNCTION COMPARF ( A , B) IS EQUIVALENT TO (A-B) BUT WILL NOT OVERFLOW. 

NORMAL INPUT - EL I ST , TBEGIN, TEND, TAPE 3 


EL I ST- THE BCD LIST OF EPHEMERIS DATA NAMES TO BE PLACED ON 

TAPE 3 . THE NAMES ARE READ FROM CARDS, AND IS USED TO 
MAKE THE T MAKE LIST. EL I ST IS NOT CHANGED IN STORAGE UNTIL 
THE FINAL CLEAR FOR THIS SUBROUTINE. 

TMAKE- THE LIST OF EPHEMERIS NAMES WITH DUPLICATES DROPPED AND 

ZERO SPACES CLOSED IN. AS THE EPHEMERIDES ARE FINISHED THE 
NAMES ARE ERRA5ED FROM THIS LIST. 

TMADE- LIKE TMAKE BUT IS HELD FOR OUTPUT. 

TBEGIN- THE BEGINNING DATE EXPRESSED AS A JULIAN DAY. 

TEND- ENDING DATE EXPRESSED AS A JULIAN DAY. 

INTVAL- THE APPROX. NUMBER OF DAYS COVERED BY ONE SET OF COEFF. IT 
IS USED TO DEC I OE WHICH DATA ARE TO BE ENTERED DOUBLE. THE 
DOUBLE ENTRIES PERMIT FASTER OPERATION IF REVERSAL OF 
INTEGRATION IS REQUIREO FOR ANY REASON. 

EDATE- JULIAN ENDING DATE FOR THE MASTER EPHEMERIS. 

ERTOAU- EARTH RADII PER A.U. 


COMMON C 






DIMENSION 

C 

( 1600), 

TMAKE 

( 12), 

LIST 

(30), 

EDATE 

( 12) , 

INTVAL 

130), 

KTAG 

(12), 

> EL I ST 

( 12) , 

TMADE 

( 12) , 

INTVA 

(2), 

PNAME 

(30) , 

TDATUM 

( 1100) , 

DATUMT 

(21,12) 
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EUUI VALENCE 

II TAPE 3 , C ! 2)1, < ERTOAU *C { 3)),( KTAG,C! A)), I FILE, Cl Ul), 

2 ( EL I ST , C I 17 ) ) » ( TBEGIN »C l 29)),! TEND, Cl 30) ) . C PNAME, C! 31J), 

31 KHAMP , C I 61)), ( THADE ,C l 73)), l TMAKE » Cl 85 ) ) , ITDATUM, C I 441 )) , 

4! EDATE ,C l 127) ) , II NTVAL , C (1 57 ) ) , ( INTVA,CU56) ) , { DATUMT , C ( 1 89 ) ) 

C 

C PART 1. REWIND 3 AND CLEAR COMMON. 

B COMPARE l A ,8) = ( A+B ) * I - I A*B > ) 

REWIND 3 
DO 1 K = L , 1600 
1 C(K) = 0.0 

C 

C THE FOLLOWING NH STATEMENTS LOAD THE BODY NAMES INTO THE MACHINE. 

C NOTE. THE EARTH IS NOT IN THIS LIST (NO EPHEMER IS FOR EARTH.) 

PNAME(l) = 3HSUN 
PNAMEI2) = 6HMERCUR 
PNAME(J) = 5HVENUS 
PNAME14) = 4HMARS 
PNAME I 5 ) * 6HJUPITE 
PNAMEI6) = 6HSATURN 
PNAME ( 7) * 6HURANUS 
PNAME I 8 ) = 6HNEPTUN 
PNAME ( 9 ) = 5HPLUT0 
PNAME (10)= 4HMD0N 
PNAME! 11) = 6HEARTHM 

C 

C PART 2. SET UP JULIAN DATES ENDING EACH EPHEMER IS . 


EDATE ( 1 ) = 2451B72.5 11/24/00 
EDATE ! 3 ) = 2451848.5 10/31/00 
EDATE ( 4 ) = 2451020.5 7/26/98 
EDATE l 5 ) * 2473520.5 2060 
EDATE ( 6 ) * 2473520.5 2060 
EDATE ( 7 ) = 2473520.5 2060 
EDATE ( 8 ) = 2473520.5 2060 
EDATE ( 9 ) * 2473520.5 2060 
EDATE (101= 2440916.5 li/26/70 
EDATE (11)= 2451848.5 10/31/00 


INTVA « 30000 
INTVAL ( 1 ) = 8 
INTVAL ( 2 ) = 5 
INTVAL ( 3 ) » 15 
INTVAL (4) * 44 
INTVAL ( 5 ) * 330 
INTVAL ! 6 ) * 825 
INTVAL! 7) * 1211 
INTVAL ( B ) * 1172 
INTVAL(9) * 1101 
INTVAL CIO) * 2 
INTVAL (11) * 15 
FILE = l. 

ERTOAU = 4.26546512 E~5 

2 END FILE 3 
MOON * 0 
LI * 1 

C 

C PART 2B. CALL INPUT AND SEE IF TAPE IS TO 8E MADE. INPUT MUST ALWAYS 

C MAKE TAPE 3 = 0. 0 IF TAPE IS TO 8E MADE. 

TAPE3 = 3. 

CALL INPUT(300»C,LIST) 

IF (TAP63) 63,3,63 
C 

C PART 3. TAPE IS TO BE MADE SO MOVE EFH6MERI S LIST TO TMAKE AND 

C TO TMADE (FOR OUTPUT), CANCEL ANY ZERO OR OUPLICATE NAMES. 

3 KOUNT = 1 
DO 6 K= 1 , 12 
TMAKE ( K ) = 0. 

TMADE ( K ) = 0. 

4 DO 5 J»l, KOUNT 

IF (COHPARF (ELIST(K) ,THAKE( J-l) )) 5,6,5 

5 CONTINUE 

TMAKE ( KOUNT ) - ELIST(K) 

TMADE ( KOUNT ) = ELIST(K) 

KOUNT = KOUNT ♦ 1 

6 CONTINUE 

KOUNT = KOUNT - l 
C 

c PART 4. FIND INPUT ERRORS. 

7 IF! TBEGIN- 2437202- 5) 66,9,9 
9 KM = 2 

11 ERROR = 0. 

WRITE TAPE 3, FILE 
DO 21 J=l, KOUNT 
KTAG(J) * 0 

12 00 13 K=1 , 20 

IF ( COHPARF (PNAHE IK) • TMAKE € J 1 1 1 13,16,13 

13 CONTINUE 
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c PART 5. PRINTS OUT THE MISSPELLED NAMES AND OTHER ERRORS. 

14 PRINT 15, TMAKE(J), TBE&IN, TEND 

WRITE OUTPUT TAPE 6 , 15, TMAKE ( J ) , TBE&IN. TEND, I PNAME ( K J , 
lEDATE(K) ,K*1,20) 

15 FORMAT ( 2 3H TROUBLE ON TAPE 3 MAKE / 2X,A6,10H T BEGIN* FlO.l.BH 
l T END* F10.1//2(2X,A6,F20.l) J 

ERROR = 1- 
GO TO 21 
C 

C PART 4B. CHECKS DATES AND STORES INDEX FOR MOON SO THAT EARTH 

C RADII CAN BE CONVERTED TO A.U. 

16 IF (10-K) 18,17,18 

17 MOON * J 

18 KTAG(J) = K 

19 IF [ EDA TE ( K ) - TEND) 14,21,21 

21 CONTINUE 

ASSIGN 36 TO NS1 
IF (ERRDR) 22,22,68 
C 

C PART 6. FIX UP A TAG ( K TAG ) TO INDICATE WHETHER TO ENTER DATA DOUBLE OR 

C NOT. KHAMP WILL BE SHORTEST INTERVAL. KT AG WILL BE NON-ZERO IF 

C ANY DATA ENTERS MORE THAN ONCE FOR 10 ENTRIES OF THE MOST 

C FREQUENT DATA. 

22 KHAMP * INTVAL(O) 

DO 23 J = l » KOUNT 
K * KTAGIJ) 

KHAMP = XMINOFIKHAMP, INTVAL (K) I 

23 CONTINUE 

KHAMP * KHAMP *10 
DO 24 J«l, KOUNT 
K * KTAGIJ) 

24 KTAG(J) = INTVAL(K) l KHAMP 
C 

C PART 7. LOCATE FILE 2 ON TAPE B. 

S 25 RTB 8 
S STZ Jl 

S CP Y DUD 

S CP Y KFILE 

S TRA *26 

S TRA *25 

S TRA *25 

26 IF (KH-KFILE) 27,32,29 

27 IF (KFILE - 3) 28,28,29 

20 CALL BKF I LE ( 8 I 
GO TO 25 

C BY PASS A FILE. 

529 RTB 8 

530 CPY DUD 

S TRA *29 

5 TRA *25 

S TRA *29 

C 

c PART 8. THIS IS CORRECT FILE ON TAPE 8, READ DATA. THERE CAN BE UP 

C TO 12 SETS OF DATA PER RECORD. A SET OF DATA IS 21 WORDS. 

31 Jl = -1 

S RTB B 

S CPY DUD 

S TRA *32 

S TRA *34 

S TRA *34 

32 Jl » Jl *1 

S CPY TDATUMU1) 

S TRA * 32 

S TRA *34 

S TRA *33 

33 Jl » Jl - l 

GO TO NS1 » < 36,46) 

34 WRITE DUTPUT TAPE 6,35, KF I LE , ( TMAKE ( K ) , K* i , KOUNT 1 

35 FORMAT U3H END OF FILE I3,67H ENCOUNTERED ON TAPE 8 BEFORE END TI 
IME SATISFIED FOR THE FOLLOWING Z1213X,A6)) 

GO TO 68 
C 

C PART 9. IS THIS A SATISFACTORY STARTING POINT, QUESTION MARK. 

C THE 1ST SET OF DATA FOR EACH PLANET MUST PRE DATE TBEGIN. 

C PART 9 IS EXECUTED ONLY ONCE. 

36 DO 42 J-LI, KOUNT 
DO 37 K*1 • Jl ,21 

IF (COMPARF (TDATUH(K) , TMAKE (J))J 37,39,37 

37 CONTINUE 
3S LI * J 

BACKSPACE 8 
BACKSPACE 8 
GO TO 31 

39 IF ( TD A TUM ( K* l ) -TDATUM ( M2) -TBEGIN) 40,40,38 

40 DO 41 KJ=1,21 

K I x K ^ K J ~ 1 

41 DATUMT (KJ, J) * TDATUM ( K 1 } 

42 CONTINUE 

IF (MOON) 43,45,43 

43 DO 44 K J*4, 21 

44 DATUMT ( K J , MOON ) * DATUMT ( KJ , MOON) *ERTOAU 

45 ASSIGN 46 TO NSl 
C 
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PART iO- PUT AWAY NEEDED OATA. TEST NAME* TIME OF BEGIN AND END. DO NOT 
WRITE TAPE 3 UNTIL TBEGIN PREDATES THE ENO OF THE FITTED 
INTERVAL. 50 REPEATS OLD DATA, 57 WRITES NEW DATA. THE NAMES 
ARE ERASED FROM TMAKE AS SOON AS THE DATA POST DATES TEND. WHEN 
ALL NAMES ARE GONE, RETURN TO INPUT 300 TO SEE IF ANOTHER 
EPHEMERIS IS TO BE CONSTRUCTED. 

46 DO 65 K- 1 , J 1 , 21 
DO 47 J 3 1 , KOUNT 

IF I COMP ARF I TOATUMt K ) ,T MAKE 1 J ) ) ) 47,48,47 

47 CONTINUE 
GO TO 65 

48 SWT = TBEGIN-TDATUMtK+lJ-TDATUHIK+2) 

IF I SWT > 49,49,52 

49 I F ( KTAG ( J ) > 50,52,50 

50 WRITE TAPE 3 , ( DATUMT ( K J , J ) , KJ=I,21) 

51 FORMAT ( IX,A6,FI0. 1) 

52 DO 53 K J* l , 2 l 
K 1 * K ♦ K J 

53 DATUMT { K J , J ) * TDATUM I K I — 1 ) 

IF U-MOON) 56,54,56 

54 DO 55 K J * 4,21 

55 DATUMT ( K J , J ) - DATUMT ( K J , J ) *ERTOAU 

56 IF (SWT) 57,57,58 

57 WRITE I APE 3 , I DATUMT I K J # J ) ,KJ=1,21) 

58 IF(TEND-QATUMTI2,J ) -DATUMT ( 3, J ) ) 59,59,65 

59 TMAKE(J) * 0 
DO 60 KK»1, KOUNT 
IF { TMAKE IKK)) 65,60,65 

60 CONTINUE 

WRITE OUTPUT TAPE 6, 61, F I LE , T6EGI N, TEND, KOUNT , ( TMADE l KK ) , 

1KK=1 , KOUNT) 

61 F0RMATI28H0EPHEHERI S COMPLETED, FILE 3 F3.,6H, FROM FL0.U3H TO 
l F10.1, 4H FOR 12, 18H BODIES AS FOLLOWS / 12(2X,A6)> 

62 FORMAT (IX,A6»7E 16.8/1 1X,7E16.8)) 

FILE = FILE ♦ l. 

GO TO 2 

63 WRITE TAPE 3, FILE 
REWIND 3 
REWIND 8 
TAPE 3 3 3. 

DO 64 J-3,1600 

64 CtJ) = 0 
RETURN 

65 CONTINUE 
GO TO 31 

66 PRINT 67, TBEGIN 
WRITE OUTPUT TAPE 6,67, TBEGIN 

67 FORMAT 1 33H TBEGIN PREDATES 2437202.5,IT IS F10.1J 
69 CONTINUE 

REWIND 8 

C END OF THE FORTRAN STATEMENTS. *••»*** 


C 

C 

c 

c 

C 


c 


c 

c 


c 


SUBROUTINE STDATA 

THIS ROUTINE CLEARS COMMON 4 TO 1300 AND LOADS A SET OF STANDARD DATA INTO 
THE MACHINE. ANY VALUES SET HERE MAY BE OVERWRITTEN BY INPUT 1 IN MAIN l. 


COMMON C 



DIMENSION 

PNAME 

( 12) , 

AMASS 

! CON 

19) , 

COEFN 

1 AK 

(3) , 

XDOT 

> REFER 

i C 

I 12) , 
I 1) 

RCRIT 


EQUIVALENCE 

1 ( STEPMX ,C t 20) ) , (CONSTU ,C ( 18) 
21 E TOL , C I 25 ) ) » ( ERL I MT , C ( 17) 
3 ( TFILE.CI 16)), I NPONG ,C I il) 
4(BGDYCD,CI8ll) ) , (MDDOUT ,C( 103) 
5 ( XDOT i C 1 161 ) ) , ( SPD , C I 2 53 ) 

6 ( OBL ATK ,C ( 39 ) ) , ( RESQRD ,C I 40) 
71 RMASS , C f 41 ) ) , ( GASF AC ,C 1 4 58) 
81 CON, C( 576) ) , ( AK , C I 2 33) 

CLEAR COMMON FROM 4 TO 1300. 

DO l J 3 4,1300 
l C(J) 3 0.0 


(30), 

NPONG 

(5), 

I 190) , 

ICC 

(4), 

I 15) , 

IND 

(3), 

130), 

AW 

(4), 



I ICC, C ( 238) ) , 

( I MODE , C ( 28)) 


( EREF.CI 

37) ) , 

( SQRDK , C ( 468) ) 


( RCRIT, CI911) ), 

( AMASS, CC881M 


( IND, C ( 791 ) ) , 

{ STEPS, C( 21)) 


l CONSUtCI 

36) ) , 

( COEFN, C(601)) 


(PNAME » C ( 82 1 ) ) , 

(REFER * C I 851 ) ) 


( OBL AT J , C ( 

38)), 

{ AW,CI261)) 


(DTOFFJ.Cf 

3in, 

I AU , C I 46 1 ) ) 
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C THE FOLLOWING NH STATEMENTS LOAD THE BODY NAMES INTO THE MACHINE. 

PNAME(l) * 3HSUN 
PNAME ! 2 ) * 6HMERCUR 
PNAME ( 3 > = 5HVENUS 
PNAME ( 4 I = 5HEARTH 
PNAME ( 5 ) * 4HMARS 
PNAME ! 6 ) * 6HJUPITE 
PNAME { 7 ) - 6HSATURN 
PNAME l 8 ) = 6HURANUS 
PNAME ( 9 I = 6HNEPTUN 
PNAME (10)= 5HPLUT0 
PNAME ( l I ) = 4HM00N 
PNAME (12)= 6HEARTHM 

C 

C FILL OUT SUN REFERENCE LIST. 

DO 2 K - 2,12 
2 REFER ( K J = PNAME 1 1 ) 

C 

C FILL OUT EARTH REFERENCE LIST. 

REFER ( 1 ) * PNAME ( A ) 

REFER! A) * 5HZER0+ 

REFER 111) * PNAME ! A ) 

C 

C LOAD THE REMAINING STANDARO DATA. 

AMI) = 0.5 

AM 2 ) = 0.5 

AK ( 3 ) = 1.0 

AMASS ( 1 ) * 1.0 

AMASS ( 2 ) = 1.0/6120000.0 

AMASS! 3 I - 1.0/406645.0 

AMASS ( 4 ) = 1.0/332488.0 

AMASS! 5 ) = 1.0/3088000.0 

AMASS ( 6 ) - 1.0/1047.39 

AMASS<7) * 1.0/3500.0 

AMASS! 8 ) = 1.0/22869.0 

AMASS! 9 ) * 1.0/18889.0 

AMASS! 10) * 1.0/400000.0 

AMASS ! 1 1 ) =AMASS!4)/81. 375 

AMASS! 12 ) * AMASS ! 4 ) ♦ AMASS ( i 1 > 

AU * 1 . A 95 Ell 
AW [ 1 ) = 1 • /6 • 

Awm*Awm*Awm 

AW!4)*AW!1) 

AW13) x i.-!AW(2)+!AW(ll+AW!4) ) ) 

BDDYCD = PNAME (4) 

COEFN! 83 ) * IE20 

CON! U = 0.2 

CON ( 2 ) « 0.2 

CDN ( 3) * 0.6 

CON 1 4 ) * 1.4 

CON ( 5) x 1.4 

CON ! 6 ) * 2.33333333 

CON ( 7) * 0.1 

CON! 8) » 0. I 

CON ( 9) * 0.5 

CONSTU = 1.0 E-6 

CONSU * IE-6 

ETOL * 0.01 

DTOFFJ * 244. E4 

EREF = IE-6 

ERLIMT * 3E- 6 

GASFAC » 20.064881 

ICC(l) = 79 

I CC ! 2 > * 79 

ICC ( 3 > » 79 

ICC (41 * 79 

IMODE * 1 

IND! 1 ) *2 

IND! 2 ) * 3 

INO! 3)*i 

MODOUT - 4 

NPONGU) » 2 

NP0NGI2) * 1 

NPDNG! 3 ) = 3 

NPONG! 5 ) = 1 

OBLATJ=l.6238E-3 

08LATK = 6.4E-6 

RCRITU) = 1.0 E+20 

RC R I T< 2 ) = 1.0 E + 8 

RCRIK3) = 6.14 E+8 

RCRIT(A) * 9.25 E+8 

RCRI T ( 5 ) = 5.78 E+8 

RCR I T! 6 ) * 4.81 E*10 

RCR I T ( 7 ) * 5.46 E*10 

RCR I T ( 8 ) = 5.17 E+10 

RCRIT19) = 8.61 E* 10 

RCR IT! 10) =3.81 E+10 

RCR I T ( 1 1 ) =1.60 E*8 
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RESQRD =4.068098877 E+13 

RMASS = i- 

SPO = 86400.0 

SQftDK - 1.32452139 E+20 

STEPMX* 100.0 

STEPS = l. 

TFILE - 1. 

XD0TC81 = 1.0 
WRITE OUTPUT TAPE 6,3 
3 FORMAT ( 7H0STDATA ) 

RETURN 

C 

C END OF THE FORTRAN STATEMENTS. 


SUBROUTINE TUDES 

THIS ROUTINE COMPUTES THE RECTANGULAR POSITION AND VELOCITY COMPONENTS 
WITH RESPECT TO THE EARTH MEAN EQUINOX AND EQUATOR OF 1950.0 FROM THE 
LATITUDE, LONGITUDE, AZIMUTH, ELEVATION, ALTITUDE, TOTAL VELOCITY, AND 
TIME. ALSO, WHEN TKICK DOES NOT EQUAL ZERO, A NON-ORAG VERTICAL STEP OF 
SIZE TKICK IS MADE IN CLOSED FORM (STATEMENTS 2 TO 4). THE INTEGRATION 
WILL THEN BEGIN AT TIME EQUAL TO TlME+TKICK WITH THE ORIENTATION SPECIFIED 
BY THE ABOVE FOUR ANGLES ANO THE COMPUTED VALUES OF ALTITUDE AND VELOCITY. 
FOR THE CLOSED FORM APPROXIMATION, A CONSTANT FLOW RATE (FLOW), VACUUM 
SPECIFIC IMPULSE (SIMP) AND ENGINE EXIT AREA (AEXIT) ARE ASSUMED KNOWN. 

THE ATMOSPHERIC PRESSURE IS TAKEN TO BE THE SEA LEVEL VALUE. 

COMMON C 
C 

DIMENSION 


1 

AMASS (30), 

ANGLES (4) 

, 

SINA 

(4), 


2 

COSA 

(4) , 

ANGLEB (4i 





EQUIVALENCE 






U 

X,C( 

45) ) , 

( Y,C( 461), ( 

z,ct 

47)), ( 

VX,C! 

42)) 

21 

VY,C( 

43) ) , 

( VZ,C< 44 ) } , t DTQFF J , C ( 

31)). t 

TOFFT * C I 

32)J 

3(ANGLES,C( 104) ) , 

( ALT , C ( 1 08 ) ) , ( 

VEL,C( 109) ), (ROTATE, Cl 459) ) 

41 

T I ME , C ( 

48) ) , 

( S I MP , C ( 5 )) » I 

RMASS, Ct 

41) ), ( 

TKICK, C( 

7) ) 

5 ( 

FLOW,C( 

33) ) , 

t STEPGD ,C ( 101) ) , (STEPNO, C( 102) ), ( 

AEXIT, Cl 

24) ) 

6 ( OBL ATN , C ( 

27) ) , 

( 8NAME,C(402) ),(RESQRD,C( 

40) ) , (OBLATJ , C t 

38)) 

7( 

AMASS, C < 0 8 1 ) ) , 

( SQRDK ,C ( 468) ) , ( 

SPD,C( 253) ) 




ALT1 = 0. 

VEL1 = VEL 
DELI = 0. 

DEL * 0. 

ASSIGN l TO NGO 

GREEN * 360. 0*(MQDF( ( DTOFFJ-24 37665. 5)/. 997269566, i. 1 + 

I MODF ( ( TOF F T*T I ME / SPO-. 71 97930 II/. 997269 566,1. > ) 

SINAfU * S1NF [ANGLES(i)/57. 29577951 

RADI US* 63 56 78 3. 28/SQRTF ( . 99 3 30 6 5 78 3 + .006 69 3421 68 5*S I NA( 1 ) “2) +ALT 
GO TO 8 

1 X * CO$A { 2 ) »COSA( l ) * RADIUS 

Y * $INA(2)*C0SA(l)‘ RADIUS 
Z * SINAI l ) “RADIUS 

IF (TKICK) 2,4,2 

2 RMASSO * RMASS 

RMASS * RMASS-FLOW*TKICK 

WRITE OUTPUT TAPE 6 , 3 , S TEPGO, STEPNO, ( ANGLES III , 1*1 , 4) ALT, TIME , VEL , 
l RMASSO, X,Y,Z 

3 FORMAT ( 6H0STEP=F5. ,2H *F4.,4X,6H L A T . = 1 PG l 5 . 8 , 7H LONG. =G15. 8, 6H AZ 
1I.*G15.8,7H ELEV.=G15.8,6H ALT . =G 1 5 . 8 /6H T I ME=G 15. 8 , 6H VEL.=GI5.8, 
67H RHASS = G15.8,4X,2HX*G15.8 , 5X , 2HY = Gl 5 . 8 , 4X , 2HZ*G1 5 . 8 I 

TIME = TIME+TKICK 

B1 = LOGF ( RMASSO /RMASS ) 

SIMPSL = SIMP-AEXIT/FLOW*10332.275 
VEL L = VEL+SIMPSL*9.80665*B1-G‘TK ICK 

ALT1 * TKICK*IVEL-G*TKICK/2.*9.80665*SIMP$L*II.-BL*RMASS/ 

1 ( RMASSO-RMASS ) ) ) 

4 RADIUS = RADIUS ♦ ALTl 

GREEN * GREEN ♦ ROTATE*TKI CK*5 7 . 2957795 
ASSIGN 5 TO NGO 
GO TO 8 

5 X * C0SA(2) »COSA ( 1 )*RADIUS 

Y = SINA(2)*C0SA(1) ‘RADIUS 
Z * SI NA ( 1 } ‘RADIUS 

IF (OBLATN-BNAME) 7,6,7 

6 DELI = ATANFCIC2-U > / < C 3- 1 . ) *S INA ( 1 ) /C3SA< 1 M *57 . 2957795- ANGLE SU) 

7 DEL2 * RADI US /G‘SINA(1)‘C0SA(1)*R0TATE‘R0TATE*57. 29 577951 
DEL * DELI + DELZ 

ASSIGN 10 TO NGO 

8 ANGLEBI1) * ANGLES I 1 ) ♦ DEL 
ANGLEB ( 2 ) = ANGLES ( 2 ) ♦ GREEN 
ANGLE613) * ANGLES t 3} 

ANGLEB ( 4 t * ANGLES ( 4) 
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DO 9 1*1,4 

SINA(I) = SINFIANGLEBII ) /57 . 2957795 ) 

9 CDSA(I) = COSF(ANGLEB(n/57. 2957795) 

Cl * 5. *R£SQR D/RAO I LiS/RAO I US*OBL AT J 
C2 * Ci*(SINA(l)*$INA(i)-*6) 

C3 = Ci*(SINA( D*SlNA(l)-.2) 

G - SQRDK* AMASS( 4 1 /RAO I US/RADIUS 
GO TO NGO, (1,5,10) 

10 C0S1 = COS A { l } *SINA (4 ) -COS A ( 4) *CQSA( 3) «S 1NA { l ) 

C0S2 = COS A (4)»SINA(3) 

VX - VELl*(CDSl*COSA(2)-COS2*SiNA(2M-y*ROTATE 
VY = VEL1*ICOS1*SINA(2)+COS2*COSA(2) )+X*ROTATE 
VZ - VELl*(5INA(l)*SINA{4}+ COS A ( l ) *COSA (3) *COSA ( 4 ) ) 
RETURN 
C 

C END OF THE FORTRAN STATEMENTS. 


C 

C 

c 

C 

c 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


SUBROUTINE ORDER 

THIS ROUTINE TAKES THE BODY LIST READ FROM CARDS AND SORTS THEM IN 

ORDER SO THAT THE DISTANCE FROM THE REFERENCE TO EACH BODY IS 

DEPENDENT UPON ALREADY COMPUTED DISTANCES ONLY. 

ELLIPSE DATA ARE READ INTO A BLOCK OF 120 STORES RESERVED FOR 

EIGHT ELLIPSES. ONE ELLIPSE IS READ INTO A 15 STORE BLOCK. 

THE SINES OF THE 3 ANGLES ARE COMPUTEO AND REPLACE THE 3 ANGLES. 

THE COSINES ARE COMPUTED AND STORED LAST If* A BLOCK. 

A BLOCK IS ARRANGED AS FOLLOWS- 

(1) * NAME OF BODY IN BCD, ONLY 6 CHARACTERS. 

(2) = NAME OF REFERENCE BODY IN BCD, SAME RESTRICTION. 

(3) = MASS OF THE BODY IN SUN HASS UNITS. 

(4) * RADUIS INSIOE OF WHICH COORDINATES WILL BE TRANSLATED TO THIS BODY. 

(5) = SEHIIATUS RECTUM IN ASTRONOMICAL UNITS. 

(6) = ECCENTRICITY OF THE ORBIT. 

(7) = SINE OF ARGUMENT OF PERIGEE. 

(8) =■ SINE OF NODES. 

19) * SINE OF INCLINATION OF THE ORBIT. 

(10) = PERIGEE PASSAGE JULIAN DAY. 

(11) = PERIGEE PASSAGE FRACTION OF DAY. 

(12) = PERIOD OF THE ELLIPSE IN MEAN SOLAR DAYS. 

(13) = COSINE OF ARGUMENT OF PERIGEE- 

( 1 4 ) = COSINE OF NODES. 

(15)= COSINE OF INCLINATION OF THE ORBIT. 

DEFINITIONS— NOTE. COMMON EXTENSION IS TRANSFERRED TO ORUM 2 DURING SEG2. 

AMASS = HASS OF EACH BODY, SUN HASSES. ORDER OF PNAME. COMMON EXTENSION. 

BMASS = SELECTED FROM AMASS. CORRESPONDS TO BNAME LIST. COMMON EXTENSION. 

BNAME = THE ORDERED LIST OF BCD BODY NAMES. CAN BE USED IN OUTPUT. COMMON. 

BODYCD = THE ORIGINAL BCD NAMES READ FROM CARDS. COMMON EXTENSION. 

BODY L = THE LIST OF BCD BOOY NAMES WITH THE REFERENCE BODY AT TOP. 

INITIALLY EQUAL TO BODY CARD LIST (BODYCD). COMMON EXTENSION. 

IBODY = ARRAY OF SUBSCRIPTS. WHEN A DISTANCE IS FOUND FROM EPHEMERIS, IT 
HAY BE ADDED (OR SUBTRACTED) FROM THE BODY POSITION GIVEN BY 
XP( IBODY) TO OBTAIN THE POSITION OF THE PRESENT BODY. COMMON. 

KZERO = COUNT OF ZERO REFERENCES. THERE MUST BE ONE AND ONLY ONE ZERO. 

NAME = ARRAY OF SUBSCRIPTS. GIVES OLD LOCATION OF NAMES IN BODYL 

FROM LOCATION IN BNAME LIST. NOT IN COMMON. 

MANE = ARRAY OF SUBSCRIPTS. INVERSE OF NAME. GIVES NEW LOCATION OF 
BNAME LIST IN TERMS DF BODYL. NOT IN COMMON. 

N8QDYS = COUNTED INTERNALY. TOTAL NUMBER OF BODYS. 

MBODYS = COMPUTED INTERNALY. TOTAL NUMBER OF EPHEMERIDES (NBODYS-1). 

NEFMRS = ARRAY OF SUBSCRIPTS. GIVES LOCATION OF BODY IN PNAME LIST 
IN TERMS OF THE EFMRS LIST. STORED IN COMMON. 

NREFER = ARRAY OF SUBSCRIPTS. LOCATES THE REFERENCE BOOY IN BODYL. 

ORDER OF THE ARRAY CORRESPONDS TO BODYL. NOT IN COMMON. 

NNREFR = ARRAY OF SUBSCRIPTS. LIKE NREFER BUT REFERS AND CORRESPONDS TO 
BNAME LIST. NOT IN COMMON. 

PNAME = A PERMANENT LIST OF BCD BODY NAMES. 1 WORD EACH <6 CHARACTERS 

MAX). USED TO IDENTIFY MASS, REFERENCE NAMES, ETC. THE LIST IS 
A MAXIMUM OF 30 NAMES. PRECISION TAPE NAMES ARE FROM 1 TO 20, 
ELLIPTIC NAMES ARE FROM 21 TO 30. COMMON EXTENSION. 

REFER = A PERMANENT LIST OF BCD BODYS THAT ARE THE REFERENCES OF 

DISTANCES GIVEN IN EPHERHERIDES (TAPES OR ELLIPSE). CORRESPONDS 
TO PNAME LIST. STORED IN COMMON EXTENSION. 

COMMON C 


DIMENSION 






1 

AMASS 

(30) , 

BMASS 

(8) , 

BNAME 

(6) , 

2 

BODYL 

(8) , 

EFMRS 

(7), 

IBODY 

(8) , 

3 

MANE 

(8) , 

NAME 

(8) , 

NEFMRS 

(8) , 

3 

NEFMRT 

(8) , 

NNREFR 

(8) , 

BODYCD 

(B), 

4 

NREFER 

(8) , 

PNAME 

(30), 

RBCRIT 

171 . 

5 

RCR IT 

(30) , 

REFER 

(30) , 

TOATA 

( 18,7), 

6 

TDEL 

(7) , 

TIM 

I7>, 

EL IPS 

( 120), 

7 

NDUD 

(9) 
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EQUIVALENCE 
It AMASS ,C ( 68 1 ) ) 
2 t BMASS ,C(417) ) 
3( BNAME » C 1 402 ) > 
41B00Y L tC ( 801 ) ) 
5 ( EFMR$,CI 410) > 
6 ( 1 BODY , C ( 425 ) I 


, ( MBOOYS »C 1 441) ) 
i (NBODYS ,C < 489) ) 
, (NEFMRS »C ( 433 ) ) 
, CBODYCO »C ( 8 1 L ) ) 
,( RMASS »C ( 41)) 
,( FILE,C(249)> 


( GK2M.CI469)) 
( GKH, C I 470 ) ) 
(P NAME, Cl 821)) 
{RBCRIT, C(450) ) 
{ RCRI T » C ( 91 1 ) ) 
( REFER* C( 85 1 ) ) 


[ SQRDK ,C 1 468 ) ) * 

( TDATA,C( 276 )) * 
l TOEL * C ( 592 ) ) , 

( TIM, C< 585)) * 

( ELIPS,C<941) ), 
(MANE ( 1 ) * NDUDI 2 ) ) 


THIS SECTION SEES WHAT ELLIPSE DATA WAS READ FROM CAROS ANO PUTS THE 
NAMES IN PLACE SO THAT DATA WILL BE USED IF NEEDED. ELLIPSE OATA HAS 
PRIORITY OVER TAPE DATA BECAUSE LAST DATA IN LIST IS THAT ACTUALLY USED. 
FUNCTION COMPARF { A , B ) IS EQUIVALENT TO (A-B) BUT WILL NOT OVERFLOW. 


C 

C 

C 


C 

C 

C 


c 

c 


c 

c 


c 


COMPARF (A, B) 
DO 3 K=l,120, 
IF(ELIPSCK) ) 

1 KOUNT = ( K- 1 ) 
PNAME I KOUNT ) 
REFERUOUNT) 
AMASStKOUNT) 
RCRI TUOUNT) 
DO 2 J = 6,8 
I*K+ J 

ELIPS(I+6} * 

2 ELIPSU) 

3 CONTINUE 


- (A+B)*(“CA*B)) 
15 

1.3,1 
/ 1 5*2 1 

- ELIPSU) 

* EUPSUd) 

- ELIPSU + 2J 

- EL I PS ( K + 3 ) 


COSFIELIPS(I) ) 
SINF(ELIPSU)) 


PART 0. THROW AWAY BLANKS AND DUPLICATES IN BNAME LIST. 

ALSO COUNT THE BODIES. 

4 DO 5 K»i,8 

5 BNAME ( K-* 1 ) * BODYCO(K) 

L * 1 

BODYL(O) * 0. 

DO B 1*1,9 
BODYL ( I ) * 0. 

DO 6 K* 1 *L 

IF (COMPARF (BNAME(I), BODYL(K-l) ) ) 6,7,6 

6 CONTINUE 

BODYL I L ) * BNAME C I ) 

L * L*l 

7 BNAME I I ) * 0. 

8 CONTINUE 
NBODYS - L-l 
MB0DY5 * NBODYS- 1 

PART 1. FIND THE REFERENCE BODY FOR EACH BODY IN THE LIST OF BODYS 
READ FROM CAROS. CLEAR NREFER ANO BNAME. 

DO 13 KL*i, NBODYS 
NREFER I KL ) * 0 
NEFMRTUL) *0 
BNAME (KL) * 0. 

DO 12 KP= 1,30 

IF (COMPARF (BODYLUll ,PNAHE(KP) ) ) 12,9,12 

9 NEFMRT(KL) = KP 

DO ii KR - 1,8 

IF ( COMPARF (REFER! KP) , BODYL IKK) ) ) 11,10,11 
10 NREFER ( KL ) * KR 
U CONTINUE 

12 CONTINUE 

13 CONTINUE 

PART 2 . COUNTS 0 REFERENCES AND SAVES TEMPORARY SET OF INDEXS. 

L4 IF (NBODYS) 24,24,15 

15 KZEROS = 0 
MISPEL * 0 

DO 20 K = 1.NB0DY5 
NNREFR(K) = NREFER ( K ) 

16 IF ( NEFMRT ( K) ) 18,17,18 

17 MISPEL = MISPEL + l 

18 IF (NREFER(K) ) 20,19,20 

19 KZEROS * KZEROS ♦ 1 

20 CONTINUE 

21 IF (KZEROS- 1) 24,22,24 

22 IF (MISPEL) 24,23,24 

23 IF (NBODYS-8) 2 8,28,24 


PART 3 . REPORTS ERRORS IN BODY LIST. 

24 WRITE OUTPUT TAPE 6,25 , NBODYS , M I SP EL , <1 EROS ,( BODYL ( K ), K* 1 , NBODYS ) 
WRITE OUTPUT TAPE 6,26 , (NREFER(K) ,K=1, NBODY S ) 

WRITE OUTPUT TAPE 6,27 , (K , PNAME ( K ), REFER ( K ) ,K= l , 30 ) 

25 FORMAT (26HOGOOFY BODY LIST (NBODYS *I2,13H, MISSPELL =12, 

1 tlH, KZER05 =12, IH) /i 1H0B0DYL 1ST =8(3X,A6)> 

26 FORMAT (11H NREFER =16,719) 

27 FORMAT ( / 5 ( 3H K3X , 4H80DY4X , 5HRE F C R5X , > /5 ( I 3 , 2X , A6 , 2X , A6 , 5X ) ) 

GO TO 50 
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C PART 4. TRACES OUT ..REFERENCE TO BODY.. RELATIONSHIPS 

28 KK =2 
KN = 1 
NAME ( l ) = 1 

29 IF (NREFERIKN)) 24,31,30 

30 NAME ( KK ) = NNREFR ( KN ) 

NNREFR { KN ) = 0 

KN * NAME ( KK ) 

KK = KK + 1 
GO TO 29 
C 

C PART 5. TRACES OUT ..BODY TO REFERENCE.. RELATIONSHIP 

31 DO 34 KN * 1 , NBODYS 
DO 34 K = 1, NBODYS 

32 IF ( NNREFR I K ) - NAME ( KN ) ) 34,33,34 

33 NAME ( KK ) = K 
KK * KK ♦ 1 

34 CONTINUE 
C 

C PART 6. INVERTS NAME TO MANE, STORES BNAHE, BMASS , RBCRIT, AND A 
C TEMPORARY NEFMRS. 

DO 35 K = 1, NBODYS 
N * NAME ( K J 
MANE (N I * K 
NEF » NEFMRT ( N) 

BNAME(K) = PNAMEINEFI 
BMASSUJ - AMASS ( NEF ) 

RBCRIT(K) « RCRIT(NEF) 

NEFMRS ( K ) * NEF 

35 CONTINUE 
C 

C PART 7. FINDS NNREFR REFERENCE FDR BNAME LIST , ALSO TEMP. I BOD Y 

DO 36 K = 1, NBODYS 
N = NAME I K ) 

NRF = NREFER(N) 

NNREFR ( K ) = MANE ( NRF } 

36 I BODY ( K J * MANE ( NRF } 

C 

C PART 8 . FINOS IBOOY FOR BACKWARD REFERENCE. 

DO 39 Mi, 8 

37 I F ( NNREFR C K ) ) 24,40,38 

38 N * NNREFR ( K ) 

IBODYIN) = -K 

39 CONTINUE 

C IBOOY LIST IS COMPLETE. 

C 

C PART 9 . WRITES OUT EPHEMER I $ LIST TO BE USED IN STORING DATA AND 
C MAKES FINAL NEFMRS LIST. 

40 KK * l 

DO 43 MI , NBODYS 

41 I F ( NNREFR ( K ) } 42,43,42 

42 EFMRS(KK) * BNAME ( K ) 

NEFMRS ( KK ) = NEFMRS ( K ) 

KK = KK + L 

43 CONTINUE 

NEFMRS ( NBODYS ) * 0 
C 

C PART 10. SAVES ELLIPSE DATA 

FILE » 0. 

DO 48 Ml, MBODYS 

44 IFINEFMRSIM-20) 47,47,45 

45 DO 46 J = 5 , 15 

L« ( NEFMRS ( K ) - 21) * 15 4-J 
TDATA I J-4, K ) - ELIPS(L) 

46 CONTINUE 
GO TO 4 B 

C 

C PART 10A. LOAOS A FALSE {VERY EARLY) TAPE TIME TO FORCE TAPE 

C READING BY THE EPHMRS ROUTINE. FILE * 0 UNLESS TAPE IS USED. 

47 TDEL(K) = 0. 

TIM(K) * 2400000.5 
FILE * 10. 

48 CONTINUE 
C 

C PART 11. COMPUTE GRAVITATIONAL CONSTANTS. 1.9866 E+30 - KILOGRAMS/SUN MASS 

GK2M = SQROKMBMASSf 1 ) +RMAS S/ 1.9866 £ + 30 I 
GKM » SQRTF ( GK2H ) 

C 

C PART 12. WRITES THE BNAME LIST ON TAPE 6 . 

WRITE OUTPUT TAPE 6 , 49 , BNAME ( 1) ,( BN AME { K ), M2, NBODYS ) 

49 FORMAT ( 19H0REFERENCE BODY IS A6,5X,23H PERTURBING BOOIES ARE 

1 7 ( 2X , A6 } } 

RETURN 

50 CONTINUE 
C 

C END OF THE FORTRAN STATEMENTS. •••••»«• 
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C MAIN 2 CONTROLS THE PROGRAM SEQUENCING FOR THE SECOND SEGMENT. IT ALSO 

C CONTAINS THE INTEGRATION SCHEMES. THE SET OF INTEGRATION VARIABLES IS 

C IDENTIFIED BY IMODE ACCORDING TO ThE FOLLOWING 
C 

C IMODE VARIABLES 

C l ORBIT ELEMENTS 

C 2 rectangular 

C 3 RECTANGULAR TEMPORARY 

C -1 ORBIT ELEMENTS— CHANGE TO RECTANGULAR 

C -2 RECTANGULAR — CHANGE TO ORBIT ELEMENTS 

C -3 ORBIT ELEMENTS— CHANGE TO TEMPORARY RECTANGULAR 

C 

COMMON C 
C 

0 1 MENS I UN 

1 XPRIM (15,2), XPRIMB (15,2), XDOTPM (15,2), 

2 X (15), XINC (15), OLOINC (15), 

3 XDOT (15), RB (3), XK (15), 

4 C ( 1 ) , AK (3) , AW ( A) , 

5 X WHOLE (15), VX (3) 

C 

equivalence 

1( IMODE, C( 28)),{TRSFER,C(22A)>,(XWH0LE,C(5AA)>,( XPRIM, C( AD), 

2 ( XPR I MB , C ( 71)), ( RATIO, C(600) ), ( XDOT , C ( 161) ), ( DELT , C ( 10)), 

3 { AH , C t 26 D ) , ( AK »C ( 2 33 ) ) , ( OL DDEL , C ( 225 ) ) , ( ACOEF 1 , C ( 265 ) ) , 

A( AC0EF2,C (266) > , ( ACOEF 3 , C ( 267 ) ) , ( X 1 NC , C ( 1 A6 ) ) , ( £2,C(260I), 

5 ( ERL IMT ,C ( 1 7 ) ) , ( KSUB , C { 2 5A) > , ( DEL , C ( 2 55 ) ) , ( STEPGO , C ( 1 0 D ) , 

6 ( ST£ PNO , C ( 102) ) ,( A SQRD * C ( 563 ) ),( GK2M, C ( A69 ) ) , ( RE VS , C ( A90 ) ) , 

7 { £ TOL , C { 25)), ( TTE ST , C 1 2 51 ) ) , ( CQNSTU, C ( 18) ) , (ASYMPT ,C( 5A3) ) , 

8 { TRU , C ( A83 ) ) , t VSQRD ,C I A76) ) , ( RB,C(200)1,( VX,C(A72)), 

9 ( ERLOG , C ( 2 59 ) ) , ( A1,C(236)),( X, C U 3D ) , ( STEPMX , C ( 20)) 

EQUIVALENCE 

U N$TART ,C (2A7 ) > , ( R , C ( AA2 ) ) , ( MBODYS , C ( AA 1) ) , ( T IME , C ( 1 38 ) ) , 

2(L£NGTH,C(257) ) , ( H2,C(256)),( A2,C(237)),( ZN,C(A87)) t 

3 ( EMONE , C ( 2 A3 ) ) 

C 

C PART 1. SET UP THE STARTING SEQUENCE FOR ERROR CONTROL AND DELAY CHECKING 

C THE ERROR UNTIL TWO STEPS ARE COMPLETED. THE ASSIGNED GO TOS NSTART AND 

C I BEG I N CONTROL STARTING. REWINDING 2 USUALLY SAVES TIME ON PING-PONG TAPE. 

REWIND 2 

1 DO 2 J = I , 8 

XPR I M( J » 2 ) = XPR I M ( J , I) 

XPR I MB ( J , 2 ) = XPRIMBUtD 

2 X( J) = XPR I M ( J , D 
NSTART * 0 

H2 - DELT 
OELT = DELT/2. 

CALL EQUATE 
CALL OUTPUT 
DO 3 J=i ,3 
XWHOLE ( J ) =VX ( J ) 

3 XWHOLEU + 3) = RB(J) 

C CHANGE INTEGRATION VARIABLES IF IMODE IS RETURN FROM TESTTR IS AT 

C BEGINNING OF MAIN 2. 

IF (IMODE) A, 5, 5 
A CALL TESTTR 

5 ASSIGN 21 TO NSTART 

C STATEMENTS 7 TO 9 INITIALIZE NREV1 AND NREV2 FOR USE IN PART 7A. 

IF ( RB ( 2 ) ) 7,6,8 

6 IF ( VX ( 2 ) ) 7,8,8 

7 ASSIGN 37 TO NREVI 
ASSIGN 35 TO NREV2 
GO TO 9 

8 ASSIGN 33 TO NREVI 
ASSIGN 37 TO NREV2 

9 DO 10 J= 1 * 8 

XDOT PM ( J , 1) = XDOT ( J ) 

XINC(J) = 0. 

10 CONTINUE 

11 KSUB = l 
ASSIGN 16 TO N 

C 

C PART 2. RUNGE-KUTTA SUBINTERVAL SCHEME. EQUATE PRODUCES THE NECCESSARY 

C DERIVATIVES XDOT(J). 

12 DO 13 J=l,8 

XK ( J ) * XOOT(J) • DELT 

XINC(J) = XINC(J) * AH ( KSUB ) *XK ( J ) 

13 XU) - XPR I M ( J , 2 ) ♦ AK ( KSUB ) *XK I J ) 

IA CALL EQUATE 

CALL DUMP (3, C, LENGTH) 

15 GO TO N, ( 16,17, 18,20) 

C 
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C PART 3. SUBINTERVALS 2, 3, AND 4* TO STATEMENT 19 FINISH A 

C RUNGE-KUTTA STEP AND INCREMENT XPRIM(J,2) IN DOUBLE PRECISION. 

16 K SUB = 2 
ASSIGN 17 TO N 
GO TO 12 

1 7 K SUB = 3 
ASSIGN 18 TO N 
GO TO 12 

18 DO 19 J = l » B 

XINCIJ) = XINCU) ♦ AW(4> *XDOT I J ) * DELT 
CALL EX ADD ( XPR I Ml J , 2 ) , XPRIMB(J,2), XINC(J>> 

XU) = XPRIMU.2) 

19 CONTINUE 

PART A. BEGIN A NEW RUNGA-KUTTA STEP. THIS ALSO GIVES DERIVATIVES 
FOR THE LOWER ORDER INTEGRATION CHECK. 

ASSIGN 20 TO N 
GO TO l A 

20 GO TO NSTART , (27*23*21) 

C 

C PART 5. STARTING PHASE PROGRAM. 

C PART 5 A . THIS SECTION COMPLETES THE FIRST STEP OF STARTING PHASE. 

21 ASSIGN 23 TO NSTART 

DO 22 1 1 8 

OLOINCI J)«XINC(J> 

XINC (J)=0. 

XDOTPM ( J » 2 ) = XDQT(J) 

22 CONTINUE 
GO TO 11 

C 

c PART 5B. MAX ERROR TE ST — ST ART I NG ONLY — CHECK THE MAX ERROR AND 

C EITHER ENTER RUNNING MODE OR REPEAT START WITH SMALLER STEP. 

23 DO 2 A J = 1 ,7 

2 A XINCIJ) * I XINC IJ ) ♦OLD INC I J ) ) *3.- 1 XDOTPM I J » 1 ) +XDQTPM IJ t 2 ) • A. 
l + XDQTU) ) *DEL T 
CALL ERRORZ 

25 IF ( E2-ERLI MT) 2 6,26,56 

26 ASSIGN 27 TO NSTART 

ASSIGN ll TO I8EGIN 

A 1 = A 2 

GO TO 31 
C 

C PART 6. RUNNING PHASE PROGRAM. 

C PART 6 A . CHECK THE INTEGRATION BY INTEGRATING OVER THE LAST 

C RUNGE KUTTA STEP BUT USE DOTS FOR LAST TWO INTERVALS, 3LDDEL 

C AND DELT RESPECTIVELY. STATEMENT 28 IS THE LOWER INTEGRATION 

C MINUS RUNGE-KUTTA INCREMENTS. ERRORZ COMPUTES THE MAXIMUM RELATIVE 

C ERROR AND STATEMENT 29 TESTS THIS ERROR AGAINST THE LIMIT VALUE. 

27 RATIO = DELT/OLODEL 
HFACT=DELT/( I. ♦RATIO) 

AC0EF1~-RATI0*RATI0»HFACT 
ACOEF2 = RATIDMDELT*3.*OLDDEU 
AC0EF3=DELT+DELT+HFACT 

DO 28 J-1,8 

28 XINCIJ) * ACOEFl*XDOTPM( J»l)+ACDEF2*XOOTPMU,2)-6.*XINC( J ) 

1 + ACOEF 3» XDOT I J ) 

CALL ERRORZ 

29 IF IE2-ERLIMT) 30,30,57 
C 

C PART 7A. LAST POINT OKAY. COUNT THE REVOLUTIONS PAST THE X-AXIS. 

C A STEP GREATER THAN 1/2 REV. MAY FAIL TO ADD IN. 

30 H2 * DELT 

31 IF ( RBI 2 ) ) 32,34,34 

32 GO TO NREV1, (37,33) 

33 ASSIGN 37 TO NREVi 
ASSIGN 35 TO NREV2 
GO TO 37 

34 GO TO NREV2, (37,35) 

35 ASSIGN 33 TO NREVI 
ASSIGN 37 TO NREV2 

36 REVS = REVS ♦ l. 

S 37 LXD IMODE , ( IMODE) 

GO TO ( 38,42,42 ) , IMODE 
C 

C PART 7B . IN ORBIT ELEMENTS. ADJUST ARGUMENT OF PERICENTER AND MEAN ANOMALY 

C TO ♦ OR - PI TO MAINTAIN ACCURACY IN SIN-COS ROUTINES. 

38 IF (EMONE) 39,42,42 

39 DO 41 J ® 3 , 6 , 3 

AD J2* I NTF ( XPR IM( J, 2) /6. 28318532 +S I GNF ( . 5 , XPR IM I J,2) >) 

IF ( ADJ2 ) 40,41,40 

40 ADJ3 * -ADJ2*6. 28125 

CALL EXADD(XPRIM(J,2) ,XPRIMB(J,2) ,ADJ3) 

ADJ3=-ADJ2*.0019353072 

CALL EXADDI XPRIMU.2) , XPR I MB I J,2> ,ADJ3) 

41 CONTINUE 
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C PART 7C • ADVANCE THE REMAINING PARAMETERS, FIND NEW STEP SUE, 

C AND TEST FOR AN ORIGIN TRANSLATION. 

42 DO A3 K= 1 , 3 

X WHOLE ( K ) *VX ( K ) 

43 XWHDLE ( K*3 ) * RBI K J 
DO 44 J = 1 1 8 

XDOT PM ( J , I ) * XDOTPM ( J , 2 ) 

XDDTPM ( J , 2 ) = XDOT(J) 

XPRIHIJ.l) = XPRI M [ J , 2 ) 

XPR I MB I J , 1 J * X PR I MB l J , 2 ) 

XINC(J) = 0. 

44 CONTINUE 
OlODEL - DELT 

45 CALL STEP 

IF (MBODYS) 46,47,46 

46 CALL TESTTR 

47 GO TO (11,11,48) , IMODE 
C 

C PART 7D . IF IN TEMPORARY RECTANGULAR COORDINATES, TEST FOR RETURN 

C TO ORBIT ELEMENTS. FIRST, E IS FOUND. IF TIME HAS NOT ADVANCED 

C SUFFICIENTLY, INTEGRATION CONTINUES IN RECTANGULAR VARIABLES (STATE. 48). 

C STATEMENT 49 DETERMINES IF KEPLERS EQUATION CAUSED IMODE * 3. IF NOT, 

C AN E CLOSE TO I CHECK IS MADE IN STATEMENT 50. IF IT DID, RECTANGULAR 

C VARIABLES WILL BE USED IF THE LIMIT IS TDO SMALL (STATEMENT 52), OR 

C IF E IS 5 OR GREATER (STATEMENT 53) OR IF THE PATH LIES CLOSE TO AN 

C ASYMPTOTE (STATEMENT 55). 

48 CALL CON VT 1 (VX,C(559>) 

E XMDDE = SQRTF ( I . +ASGRD/GK2M* { VSGRD/GK2M-2./R ) ) 

EMON E* EXMODE- 1 . 

IF ( (TIME-TTEST)*DELT ) 11,11,49 

49 IF (ASYMPT) 51,50,51 

50 IF < ETOL-ABSFI EMONE ) ) 55,U,I1 

51 IF(EMONE) 55,55,52 

52 IF(C0NSTU-l.E-7) 11,53,53 

53 IF (EXMODE-5.) 54,11,11 

54 CALL C0NVT2 

IF (ABSF (TRUJ-2.2/SQRTF (EXMODE J ) 55,55,11 

55 ASYMPT = 0.0 
IM0DE--2 

GO TO 46 
C 

C PART 8. COMES HERE WHEN ERROR TEST F AILED--BOTH STARTING AND RUN. 

C RETRIEVE OLD POINT AND RECOMPUTE WITH SMALLER INTERVAL. 

C IF TWO CONSECUTIVE TRYS FAIL (STATEMENT 59) THE STARTING SEQUENCE OCCURS. 

56 ASSIGN 1 TO IBEGIN 

57 DO 58 J* l , 8 

XPR I M( J , 2 ) * X PR I M < J » l ) 

XPR I MB ( J , 2 ) » XPRIMBUfl) 

XDOT (J)-XD0TPM(J,2) 

XINCIJ)* 0. 

58 CONTINUE 

STEP NO* STEPNO+I. 

H2 - DELT 

DELTAS I GNF ( EXPF ( < E RL 00- A2 ) / 5 . ) , DELT ) 

A2 * A l 

59 IF ( FA I L-STEPGO J 60,61,60 

60 FAIL * STEPGO 

GO TO IBEGIN, (11,1) 

61 ASSIGN l TO IBEGIN 

IF (STEPNO ♦ STEPGO - STEPMX) 62,62,45 

62 GO TO IBEGIN, (11,1) 

C 

C END OF THE FORTRAN STATEMENTS. ••••»••• 


SUBROUTINE EQUATE 

THIS SUBROUTINE IS CALLED FROM MAIN 2 TO EVALUATE THE DERIVATIVES OF THE 
VARIABLES OF INTEGRATION. EITHER RECTANGULAR COORDINATES DR ORBIT ELE- 
MENTS MAY BE USED AS THE VARIABLES OF INTEGRATION, BUT IN THE CASE OF THE 
LATTER, THE CORRESPONDING RECTANGULAR COORDINATES MUST FIRST BE FOUND. 

THIS IS DONE AT THE BEGINNING THRU THE USE OF KEPLERS EQUATION. THE 
PERTURB AT I NG ACCELERATIONS ARE FOUND BY CALLING VARIOUS OTHER SUBROUTINES 
AND THEIR SUM RESOLVED ALONG THE X,Y,Z AXIS. FINALLY, THE DERIVATIVES 
ARE CALCULATED. IN THE CASE OF ORBIT ELEMENTS, THE X,Y,Z PERTURBAT ING 
ACCELERATION COMPONENTS MUST FIRST BE RESOLVED INTO CIRCUMFERENTIAL, RADIAL 
AND NORMAL COMPONENTS. THIS ROUTINE ALSO CHANGES THE INTEGRATION VARI- 
ABLES FROM ORBIT ELEMENTS TO RECTANGULAR VARIABLES IF THE ECCENTRICITY 
APPROACHES UNITY. 

COMMON C 
C 


1 

DIMENSION 

C 

(1), 

VX 

131, 

QX (3), 

2 

RB 

(3) , 

NEFMRS 

(8) , 

X (3), 

3 

XPR I MB 

(15,2), 

FORCE 

(3), 

X IFT (3), 

4 

DRAG 

(3), 

OBLAT 

(3), 

COMPA (3), 

5 

XDD 

(6) , 

XDOTTR 

(6) , 

XPRIHC 15,2 ) 


C 
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c 

s 

c 

c 

c 

c 


c 

c 

c 


c 

c 


c 

c 


c 

c 

c 


EQUIVALENCE 
i( DM,C(L6l)J 
2 ( ASQRD f C ( 563 ) ) 
3(NSTART,C(247) ) 
4( C I NCL » C { 495 ) ) 
5 ( C I RCUH , C ( 54 L ) ) 
6 ( 5 IMP ,C ( 5)) 

7 ( CDMPA ,C ( 537 ) ) 
8 ( BNAME ,C ( 402 > I 
9(Z0RMAL,C(542>> 
EQUIVALENCE 
1 ( ASYHPT,C( 543) ) 
2 ( CONSTU , C ( 1 8 > ) 
3(CQSTRU»C(493) ) 
4 ( COSV * C ( 497 ) ) 
5 ( DE ,C { 1 62 ) ) 

6 ( ZN , C ( 487 ) ) 

7 ( D I NCL ,C ( 165) ) 
81 DNODE,C( 164) ) 
9(D0MEGA,C( 163) > 


( DMA ,C ( l 66 ) ) 

< E , C ( 1 32 ) ) 
t PRESS, C(466>) 
( EMONE ,C (243) ) 
( EPAR »C l 245 ) I 
( ETOL.CI 25>) 
(EXMODE i C ( 244 1 ) 

< FORCE ,C(525) » 
( GKM ,C ( 470) ) 

( GK2M , C ( 469 ) I 
( I MODE , C ( 28)1 

< KSU8 »C ( 2 54 ) ) 
(SINTRU, C(492) > 
( MBODYS i C ( 44 1 ) ) 
( OMEGA ,C ( 1 33) ) 
( NEFHRS , C ( 433 ) ) 
( OBL AT , C ( 5 34 ) ) 
( OBLAT J *C ( 38)) 


( P,C(137)), 
( PH I , C( 485 ) ) , 
( QX , C ( 522 ) ) , 
< R, C ( 442 ) ) , 
(RADIAL, C(540)), 
( OBLATN, C( 27)1, 
( RB, C ( 200 ) ) , 
( TOFFT, C ( 32)), 
I RMASS , C( 1 3 L ) ) , 

{ RSQRD, C( 567) ) , 
( S INCL, C ( 494) ) , 
( SINV,C( 496) ) , 
( SPD, C ( 253 ) ) , 
( DP, C ( 167 ) ) , 
( TA6LT , C( 252) ) , 
{ PUSH, C ( 34)), 
( FLOW, C( 33)), 
( TIME, C( 138) ) , 


( DRAG , C ( 531 ) ) , 
( TRSFER , C ( 224 ) ) , 
( RATMOS , C ( 248 ) ) , 
( TTEST , C ( 2 51 ) ) » 
( ZN0DE,C( 134) ) , 
( V , C ( 475 ) ) , 
( VSQRD , C ( 476 ) ) , 
( VX , C ( 472 ) ) , 
I X , C( 1 35) ) 

( XDD ,C(162)), 
( XDOTTR , C ( 132) ) , 
{ XIFT , C { 528 ) ) , 
( XPR I M , C ( 41)1, 
( XPR IMS * C ( 71)), 
( XWHOLE , C ( 544) ) , 
( ZINCL,C( 135)1 , 
( ZM, C { l 36) ) , 
( AEX I T , C ( 24)) 


TABLT*TIM£/SPD+TOFFT 
LXD IMODE , ( I MODE ) 

1 GO TO (2,16,16) , IMODE 


STATEMENTS 2 TO 16 FIND THE RECTANGULAR POSITION AND VELOCITY FROM ORBIT 
ELEMENTS AND TRUE ANOMALY. THE TRUE ANOMALY IS FOUND FROM ITERATIVE 
SOLUTION OF KEPLERS EQUATION. 

2 E2 * E*E 
E2M l x l. -E2 
EMONE x E- 1 . 

EPAR=SQRTF(ABSF(E2Ml) ) 

VCIRCL*GKM/ SQRTF ( P ) 


COMPUTE SINE AND COSINE OF TRUE ANOMALY. 
PART A. E* 1 

3 IF (EMONE) 10,4,5 

4 SINTRU = 0. 

COSTRU = 1. 

GO TO 14 


PART B. E IS GREATER THAN 1 

5 DO 7 J*l,lOO 
DELM*ZM-U*E*SINHFCU) 
ECOSU-E*COSHF(U) 

DELU * DELM/( 1.0-ECOSU) 

U ■ U+DELU 

6 IF ( ABSF(DELM)-CONSTU) 9,9,7 

7 CONTINUE 
ASYMPT * 1.0 

IF (MBODYS) 8,23,8 

8 CALL EPHMRS 
GO TO 23 

9 C05U * COSHF(U) 

DEMI » 1 • 0-E*COSU 
COSTRU - ( COSU-E) /DEMI 
SINTRU *-EPAR*SINHF(U)/OEMl 
GO TO 14 


PART C. E IS LESS THAN l 

10 DO 12 J* 1 , 5 
DELH x ZM-U+E*SINF(U) 

ECOSU « E*COSF ( U) 

DELU » DELM/( 1.0-EC0SU*0.0i*ECDSU**3) 
U = U+DELU 

11 IF 1 ABSF(DELM)-CONSTU) 13,13,12 

12 CONTINUE 

13 COSU = COSF(U) 

DEMI = 1.0-E*COSU 
COSTRU » ( COSU-E ) /DEMI 
SINTRU = EPAR#SINF(U)/DEM1 

14 PD VR * 1 • +E *COS TRU 


COMPUTE POSITION AND VELOCITY FROM ORBIT ELEMENTS AND TRUE ANOMALY. 
ALSO, CLEAR THE PERTURBATING ACCELERATIONS. 

15 SQMEGA=SINF { OMEGA) 

COMEGA*COSF (OMEGA) - 
SNODE=SINF( ZNODE ) 

CNODE=COSF (ZNODE) 

SINCL-SINFUINCL) 

CINCL*COSF(ZINCL) 

SINV=SINTRU*C0HEGA*C05TRU*S0MEGA 
COS V*C OS TRU*C OMEGA - S I NTRU*S OMEGA 
AR*COS V«CN00E-S I NV *S NOD E*C INCL 
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B 1 = S IN V*CNODE + COS V *SNQDE*C I NCL 
C 1*C OS V* SNQDE + S I NV *CNODE*C I NCL 
Dl a *SINV*SNODE-COSV*CNODE*CI NCL 
E1=E*S0MEGA+SI NV 
F i-E *COMEGA+CQSV 
AS*E i*CNODE+F l*SNODE*CINCL 
B2*F 1*CN0DE*C INCL-E1*SN0DE 
R * P/PDVR 
RSQRD = R*R 
SINVY*$INV*SINCL 
R B ( 1 ) = R* AR 

RB ( 2 ) = R*C1 

RB(3) - R *S I NVY 

VX( l }*-VCIRCL*AS 
VX(2 >=VCIRCL*B2 
VX(3)=VCIRCL*F1*SINCL 
GO TO 18 
C 

16 DO 17 K= 1 , 3 

VX { K ) =XDO TTR ( K ) 

17 RB ( K ) = X(K) 

RSORD = RB Cl) *RB ( 1 ) + RB(2)*RB(2) + RB ( 3 ) *RB ( 3 ) 

R=SQRTF (RSQRD) 

18 VSQRD=VX( l)*VX( U + VXIZ) * VX f 2 ) + VX ( 3 ) *VX ( 3 ) 

V * SQRTF ( VSQRD > 

DO 19 1=1,15 

19 C( 1+521) * O* 

C 

C TEST FOR PRESENCE OF PERTURBING BODIES. 

IF (MBQDYS) 20,21,20 

20 CALL EPHMRS 

21 IF (XABSF(IMOOE)-l) 26,22,26 

TEST FOR CHANGE FROM ORBIT ELEMENTS TO TEMPORARY RECTANGULAR 
COORDINATES IF E IS TOO NEAR TO UNITY. 

22 IF ( ETOL-AB SF { EMONE ) ) 26,23,23 

23 IF IIMODE) 54,24,2* 

24 I MODE = - 3 
IF ( NSTART ) 25,54,25 

25 TTE S T = T I ME 
CALL TESTTR 

TEST FOR OBLATENESS PERTURBATION COMPUTATION. 

26 CLA OBLATN 
CAS BNAME 
TRA • 30 
TRA *29 
TRA *30 

29 CALL OBLATE 

TEST FOR PRESENCE OF THRUST. COMPUTE THRUST MAGNITUDE IF NOT SPECIFIED. 

30 DM = -FLOW 
IF (R-RATHOS) 31,31,32 

31 CALL ICAO 
GO TO 33 

32 PRE SS = 0 . 

33 IF ( SIMP ) 34,35,34 

34 PUSH = SIMP*FL0W*9.B0665 - AEX I T*PRESS * 1 00 . 

35 I F l PUSH ) 37,36,37 

36 ASSIGN 40 TO NDONE 
GO TO 38 

37 CALL THRUST 
ASSIGN 41 TO NDONE 

TEST FOR EXISTENCE OF ATMOSPHERE. FIND AERODYNAMIC FORCES. 

38 IF (PRESS ) 39,42,39 

39 GO TO NDONE, 140,41) 

40 CALL THRUST 

41 CALL AERO 

C SUM COMPONENTS OF THE PERTURBING ACCELERATION. 

42 DO 43 J- 1 , 3 

43 COMPA ( J ) = -QX< J)+OBLATI J ) +FORCE ( J ) +X I FT( J ) +DRAGU > 

44 GO TO (47,45,45), I MODE 
C 

C COMPUTE DERIVATIVES FOR THE RECTANGULAR VARIABLES OF INTEGRATION. 

45 DO 46 K= 1 , 3 

XOD(K) * C0MPAIK)-GK2M*X(K) /R/RSQRD 

46 XDD ( K+3 ) = XDOTTRIK) 

GO TO 54 
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C COMPUTE THE DERIVATIVES OF THE ORBIT ELEMENTS. (AFTER RESOLVING 

C PERTURBATING ACCELERATION INTO CIRCUMFERENTIAL, RADIAL, NORMAL COMPONENTS) 

47 C IRCUM=COMPAI 3)*COSV*SINCL-COMPA( I > *B l -COMPA ( 2 ) *D1 
RAOIAL X COMPA(1>«AR + COMPA(2) •CUCOMPAI 3)*S1NVY 

ZORMAL = COMPAI i) »SNODE*S INCL-COMPAI2 )*CNODE*S INCL+COMPA! 3) »C INCL 

ZN*VCIRCL*E2M1»EPAR/P 

RDVPP1 = l./PDVR ♦ 1. 

RDVA = E2M1/PDVR 
DP=2.*R/VCIRCL«CIRCUM 
I F { E ) 4 B , 48 , 49 

48 CSQRD * C I RCUM*C I RCUH 
RASQRD * RADIAL*RADI AL 

DEMI = l 4. * CSQRD* RASQRD ) *VC I RCL 
VDV2R x VCIRCL/R/2* 

DE = SQRTF(4.*CSQRD+RASQRD) /VCIRCl 
DOMEGA * VDV2Rf«2.*CSQRD+RAS0RD)/DEMI*RADIAL 
DMA = ZN-VOV2R-M6. *CSQRD + R A SQRD ) /DEM 1 *RAD I AL 
GO TO 50 

49 DE * (SINTRU»RADIAL+!PDVR-RDVA)/E*CIRCUM)/VCIRCL 

DOME GA= ( SINTRU/E*RDVPP1»CIRCUM-C0STRU*RADIAL/E ] /VC I RCL 
DMA=ZN+E PAR/ VC IRCL* t < COS TRU/E-2 . / PDVR ) »RAD I AL- ! S INTRU/E*RDVPP 1»C IR 
i CUM)) 

50 IF ( S INCL ) 51,52,51 

51 DNODE = S I NV/S I NCL *ZORMA L/ VC IRCL/PDVR 
GO TO 53 

52 DNODE = 0.0 

53 DINCL « COSV*ZORMAL/PDVR/VCIRCL 

54 RETURN 
C 

C END OF THE FORTRAN STATEMENTS. **«***•• 


SUBROUTINE ERRORZ 
C 

C THIS SUBROUTINE COMPUTES THE RELATIVE ERRORS BETWEEN THE R-K AND LOW-ORDER 

C INTEGRATION SCHEMES. IT ALSO COMPUTES THE ERROR COEFFICIENT, A, AND SAVES 

C THE ERROR DATA WHEN EREF HAS A - SIGN. THE BRANCH ON IMODE DETERMINES 

C WHICH SET OF NORMALIZING FACTORS ARE TO BE USED. 

C 

COMMON C 
C 

DIMENSION RELERR ( 7 ) 

C 

EQUIVALENCE 

LI RM ASS , C ( 56)), ( E,C( 57)), ( AS, C ( 1 5 1 > ) , I OMEGAS , C ( 14B ) ) , 

2 ( RMA SSS » C ( 146) ) ,( P,C( 62)), ( ES , C < 147 > ] , I ZNODES , C ( 149 J ) , 

31 R » C ! 442 ) ) , I PS , C ( 1 52 ) ) , ( Z I NCL S , C t 150) ) , IXINC ,C(146)>, 

4 ( V»C(475) ) * I IMODE, C( 28)),( T I ME , C I L 38 H , { E2,C(260)>, 

5 ( VX , C ( 1 47 ) ),( VY , C ( 148) )•( VZ,C(149)J,( X,C(150)), 

61 Y,C(151)),( Z ,C( 152) ) , I RELERR, C< 146) ) , I A2,CI237)), 

7 ( DEL T » C ( 10)), I A 1 , C < 2 36 ) ) , { EREF , C< 37 > J , I STEPGD , C I 1 0 IJ) , 

8 ( ST£ PNO , C ( 1 02 ) ) , I I NOERR , C ( 49 1 ) ) 

C 

E2 * 0. 

RELERR! 1)*RMASSS/RMASS 
IF ! I MODE- 1 ) 2,1,2 
C 

C COMPUTE THE NORMALIZED INTEGRATION ERRORS FOR THE ORBIT ELEMENTS. 

1 RELERR I 2 )=ES/ l E+l. 0 ) / 10 • 0 
RELERR! 3) =OMEGA S/62. 831 853 
RELERR(4)=Z NODES/ 62. 831853 
RELERR! 5 )=Z I NCL S/62. 831 853 
RELERR ( 6 ) * AS/62.831853 
RELERR! 7)*P S/P/10. 0 

GO TO 3 
C 

C COMPUTE THE NORMALIZED INTEGRATION ERRORS IN RECTANGULAR VARIABLES. 

2 VI * V+100. 

RELERR! 2>=VX/Vl 
RELERR! 3)=VY/Vi 
RELERR(4)=VZ/V1 
RELERR(5)=X/R 
RELERR !6 )-Y/R 
RELERR! 7)=Z/R 

C 

C SELECT MAXIMUM ERROR, COMPUTE ERROR COEFFICIENT, POSSIBLY SAVE ERROR DATA. 

3 DO 5 J-1,7 

IF (ABSF (RELERR! J> )-E2 > 5,5,4 

4 K = J 

E 2 * ABSF ! RELERR! J) ) 

5 CONTINUE 

E2 * E2 ♦ 2E-B 
Al * A2 

A2 * L0GF!E2)-5.*L0GF( ABSF(DELT) ) 

IF (EREF) 6,7,7 

6 WRITE TAPE 4 , K , RELERR , E2 , A2 , DE L T, T I ME, STEPNO, STE PGO 
INDERR = INDERR ♦ 1 

7 RETURN 
C 

C END OF THE FORTRAN STATEMENTS. ••*•**** 
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SUBROUTINE STEP 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


c 

c 


s 


c 

c 


SUBROUTINE STEP TESTS TOR THE END Of THE PROBLEM, COMPUTES STEP SIZE, AND 
CONTROLS QUANTITY OF OUTPUT DATA. WHEN END OF PROBLEM IS DETECTED, OUTPUT 
OCCURS, THE ERROR DATA TAPE IS REWOUND, AND THE FIRST SEGHENT IS CALLED TO 
ALLOW INPUT. FOLLOWING IS AN EXPLANATION OF CONTROL ON QUANITY OF OUTPUT. 


MODOUT = I OUTPUT EVERY NTH STE P { N=STEPS ) UNTIL TIME - THIN, THEN 
GO TO MODE 2 . 

2 OUTPUT AT INTERVALS OF OELMAX UNTIL TIME * TMAX* 

3 OUTPUT AT INTERVALS OF DELMAX UNTIL TIME » TM IN, THEN 

GO TO MODE 4 . 

A OUTPUT EVERY NTH STEP UNTIL TIME * TMAX. 

COMMON C 

DIMENSION NPONGI5) 


EQUIVALENCE 

II DELT.CI 10)1,1 E2 »C 1 260) ) , I 
21 DEL ,C{ 255) ) • { ERLOG ,C < 2 59) ) , { 
3 (DEL MAX, Cl 23 ) ) , I STEPNO , C ( 1 02 ) ) , ( 
4( STEPMX ,C I 20) ) , I STEPGO , C ( 1 01 ) ) , I 
5(M0D0UT,C( 103) ) , I A2,C(237>),( 


NPONG, C I ID) 
TIME, Cl 136) i 
STEPS, Cl 2D) 
TMAX , C C 30 ) ) 
RATIO, CI600D 


A 1 , C I 2 36 ) ) , 
T M I N , C { 22)), 
SPACES, CI258D, 
H2,C(256) ) , 
TT0L.CI226)) 


PART 1. TEST FOR END OF THE PROBLEM (MAXIMUM PROBLEM TIME OR MAXIMUM 
NUMBER OF STEPS). 

STEPGO » STEPGO * l. 

IF ( ABSF ( TKAX-TIME J -TTOL ) 1,1,3 

1 CALL OUTPUT 

WRITE OUTPUT TAPE 6,2 
PRINT 2 

2 FORMAT ( 25H OCAS E COMPLET ED, T IME« TMAX 3 
GO TO 6 

3 IF < STEPGO^STEPNO-STEPMX) 7,9,* 

A CALL OUTPUT 

WRITE OUTPUT TAPE 6,5, STEPMX 

5 FORMAT (22H0ST6PG0*STEPNQ*STEPHX*F6. ) 

6 REWIND * 

CALL PONGI NPONG ( 5 ) ) 

PART 2. COMPUTE STEP SIZE ( D6LT ) AND CONTROL OUTPUT. 

7 A3 - ( A2-AL)*RATIO»A2 

8 DELT * SIGNF(EXPF( (ERLOG- A33/5. ) , DELT ) 

IF ( DEL T/H2-3 . ) 10,10,9 

9 DELT * 3 * *H2 

10 LXD MODOUT,! MODOUT) 

GO TO ( 11, 15, 13, 2i) , MODOUT 

11 IF(DELT*(TIME ♦ 3.*0ELT-TMIN)I 21,12,12 

12 MODOUT * 2 

DEL * TM IN - TIME 
GO TO 16 

13 IF ( DELT • (TIME - THIN) ) 15,15, 14 

14 MODOUT * 4 
GO TO 21 

15 DEL » DEL-H2 

16 SPACES * INTFHDEl/DELTDSIGNFl.9,(DEL/DElT) )) 

17 IF ( SPACE S ) 20, 18,20 

18 CALL OUTPUT 
DEL * DELMAX 

IF ( ABSF ( DEL ) - ABSF { DELT ) ) 19,16,16 

19 OELT = SIGNF ( DEL, DELT ) 

GO TO 16 

20 DELT * DEL/SPACES 
GO TO 23 

21 IF (MODFCSTEPGO, STEPS)) 23,22,23 

22 CALL OUTPUT 

23 GO TO (26,24,26,24) , MODOUT 

24 IF l ( TIME ♦ DELT - TMAX) *DELT) 26,25,25 

25 DELT = T MAX-T I ME 

26 RETURN 

END OF THE FORTRAN STATEMENTS. **** 


SUBROUTINE TESTTR 
C 

C SUBROUTINE TESTTR MAY BE CALLED FOR ONE OF TWO REASONS, (1) TO TEST FOR AND 
C POSSIBLY TRANSLATE THE ORIGIN (WHEN I MODE IS ♦) OR (2) TO CHANGE THE 
C VARIABLES OF INTEGRATION (WHEN IMODE IS -). A TRANSLATION OF THE ORIGIN 
C OCCURS WHEN THE OBJECT MOVES INTO A SPHERE OF INFLUENCE WHICH IS SMALLER 

C THAN ANY OTHERS IT MAY ALSO BE IN. WHEN THIS HAPPENS, THE NAME OF THE NEW 

C ORIGIN IS MOVED TO THE BEGINNING OF THE BNAHE LIST AND THE FIRST SEGMENT 

C CALLED TO REORDER THE BNAHE LIST. 

C 

COMMON C 
C 
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DIMENSION BMASSI 8 ) , BNAME (8) • R8(3,8), RBCRIT { 8 } * RREU8), CC1) 
1X13) ,XPRIM< 15*2) » XPR I H8 ( 15,2) » X WHOLE (6),VEFH(3,8), NPONGI 5 ) , 

2VX( 3 ) , OR BEL S (6) 


C 


C 

s 

c 

c 


c 

c 


c 

c 


c 

c 


EQUIVALENCE 

1C BMASS * C 1 417 ) ) i ( BNAME ,C < 402 )), l CHAMP, C< 246)) , < NP0NG,C( 111), 
2( GK2M, C ( 469 ) ) , t IMODE, C( 28) J , ( N80DYS, C ( 489) ) , < STEPNO , Cl 257 )) , 
3 ( RB,C(200)1, (RBCRIT, C(450)),( RREL » C ( 442 ) ) , ( SQRDK , C ( 468 ) > , 
4( QRBELS»C ( 227) ) , ( TRSFER ,C < 224) ) , ( X,C(200>),( XPRIM.CI 41) ), 
5 ( XPR I MB , C ( 71) ) , I X WHOLE »C ( 544) ) , ( TTEST . C ( 25 1) ) , ( VEFH,C(498)) f 
6( VX«C(472})p( REVS ,C ( 490) ) , { DELT , C ( 10)J,< TMAX,C( 30)), 


7 ( TIME ,C( 138 )). ( 


E, Cl 227) ) , ( TRU,C(483)), ( ASYMPT ♦ C ( 543 ) ) 


LXD IMODE, ( I MODE) 
IF (IMODE) 12,12,1 


IF IMODE IS TEST FOR TRANSLATluN OF THE ORIGIN, 

1 ASSIGN 27 TO N 
CHAMP* l.EOO 

DO 4 J&= 1 » NBOOYS 

IF (RREL( JB)-RBCRIT(JB) ) 2,4,4 

2 IF ( CH AMP-RBCRI T ( JB ) ) 4,4,3 

3 CHAMP * RBCRIT(JB) 

NCHAMP * JB 

4 CONTINUE 

IF ( NCHAMP- l ) 26,26,5 

5 TRSFER = 1.0 
ASSIGN 29 TO N 

8 BTEMP = BNAME (1) 

BNAME ( l ) = BNAME (NCHAMP ) 

BNAME ( NCHAMP ) * BTEMP 
TTEST * 0. 

REVS = 0. 

9 PRINT 10, BNAME ( NCHAMP ) , BNAME ( 1 ) 

WRITE OUTPUT TAPE 6 , 10 , BNAME {NCHAMP ), BNAME(l) 

10 FORMAT (28H00RIGIN IS TRANSLATING FROM A6,4H TO A6 ) 
CALL EPHMRS 

DO il K * l , 3 

VX ( K ) = VX{K)-VEFM(K, NCHAMP) 

X ( K ) = RB ( K , NCHAMP ) 

XPRIM(K+1,1)*VX(K) 

XPRIM(K*4, 1 ) * X ( K ) 

XPRIHBIK^ltl) * 0. 

XPR I MB ( K + 4 , 1 ) * 0. 

X WHOLE ( K ) * VX t K ) 

11 XWHOLEU + 3) * X { K ) 

GO TO 20 


IF IMODE IS CHANGE THE VARIABLES OF INTEGRATION. 

12 ASSIGN 28 TO N 
DO 13 K*i,3 

XPRIM{K + l,l) * X WHOLE ( K ) 

XPRIM(K>4,1) *XWHOL E { K + 3 } 

XPRIMB1 K*l, 1 ) * 0. 

XPRIMB(K+4, 1) - 0. 

VX { K ) = XWHOLEfK) 

13 X ( K ) * X WHOLE ( K+3 ) 

GO TO { 16,14,15) , IMODE 

14 CODE * 5H0RBIT 
IMODE * l 

GO TO 18 

15 IMODE * 3 
GO TO 17 

16 IMODE * 2 

17 CODE * 6HRECTAN 
ia NCHAMP * l 

PRINT 19, CODE 

WRITE OUTPUT TAPE 6, 19, CODE 

19 FORMAT ( 33H0I NTEGRATIQN MODE IS CHANGING TO A6) 

20 GO TO (21,26,26), IMODE 

21 CALL CDNVT1 (VX,CI559) ) 

GK2H* SQRDK* ( BMASS < NCHA HP ) + XPR I M l 1, l)/l. 9866 E*30) 

CALL C0NVT2 

IF ORIGIN TRANSLATION CAUSES PATH TO LIE NEAR AN ASYMPTOTE, CHANGE 
INTEGRATION VARIABLES TO RECTANGULAR IF THEY ARE ORBIT ELEMENTS. 

IF (E-l.) 24,24,22 

22 IF (ABSF (TRUI-2.3/SQRTF (El ) 24,24,23 

23 ASYMPT * 1.0 
GO TO 15 

24 DO 25 J*l,6 

25 XPRIMU*i,l) = ORBELS(J) 

26 GO TO N, (27,28,29) 

27 RETURN 

28 CALL PONG (NPONG( l ) ) 

29 CALL PONG (NP0NGI5) ) 


END OF THE FORTRAN STATEMENTS. 
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SUBROUTINE ICAO 


SUBROUTINE ICAO DETERMINE S THE ATMOSPHERIC TEMPERATURE , PRESSURE, AND 
DENSITY AS A FUNCTION OF ALTITUDE ABOVE AN OBLATE EARTH IN ACCORDANCE WITH 
NACA 1235 AND U.S. EXTENSION TO THE ICAO STANDARD ATMOSPHERE. A SHORT SAP 
PROGRAM FOLLOWS ICAO WHICH PROVIDES A MEANS OF LOADING OATA INTO MACHINE* 
IT MUST BE LOADED DIRECTLY AFTER ICAO. IF THE LENGTH OF ICAO IS CHANGED, 
THE DATA MUST 8E RELOCATED. 


R IS DISTANCE TO CENTER DF EARTH IN METERS. 

ALT IS VEHICLE ALTITUDE ABOVE AN ELLIPTIC EARTH IN METERS. 

GEO H IS THE GRAVITATIONAL POTENTIAL IN METERS. 

TABLE H IS METERS OF ALTITUDE FROM THE EARTHS SURFACE AND IS 
THE ARGUMENT OF ATMOSPHERE PROPERTY TABLE. 

ALM IS THE MEAN SLOPE DF THE TABLE H VS. TM CURVE AT TABLE H. 
TMR IS TM AT TABLE H. 

REF P IS THE PRESSURE IN MILLIBARS AT TABLE H. 

TM IS THE TEMPERATURE TIMES STD. MOLECULAR WEIGHT / ACTUAL 
MOLECULAR WEIGHT. DEGREES KELVIN. 

PRESS IS PRESSURE IN MILLIBARS. 

DNS I T Y IS DENSl'TY IN KILOGRAMS PER CUBIC METER. 


COMMON C 

DIMENSION TABLE HI 11 ) * TMRI lit, REF P(ll),ALHUl) 


21 GEO H »C t A65 I ) , I PRESS , C (466 1 ) , I TM,Ct A671 ) » (ONSI TY ,C( 4601 ) * 

3 ( TABLT,Ct252) ) , < ALT ,CI A63) I • I R,C(AA2)),( Z f C(137M» 

A I TABLE H ( 1 2 I , TMR) , (TABLE H ( 2 3) , ALM) , { T ABLE HI34),REF P) 

' ALT - R-63567B3. 28/ SORT F I . 9933065783+ . 0Q6693421685(Z/R)**2) 

GEO H = ALT/ ( l • 0 + ALT/6356766.0) 

C 

C FIND THE GEOPOTENTIAL HEIGHT IN A TABLE OF BASE DATA. DATA ARE 

C ARRANGED IN DECENDING GEO H WITH TEN REGIONS, AN UTH IS GIVEN 

C FOR EXTRAPOLATION. ABOVE THAT, PRESSURE AND DENSITY ARE SET =0. 

S LXO K, (K) 

1 IF ( K-l i } 2,6,6 

2 IF (GEO H - TABLE H ( K+ 1 ) ) 5,3,3 

3 K - K*l 
GO TO I 

A K = K-l 

5 IF <K) 7,7,6 

6 H INC = GEO H -TABLE H ( K ] 

IF (H INC) A , B , B 

7 K * I 

8 GO TO 19, II, 9,11,9,11, 9, 9, 9, 9,12), K 
C 

C CONTROL COMES HERE FOR NON I SOTHE RH AL LAYERS 

9 TM * TMR(K) ♦ ALM { K ) »H INC 

PRESS* REF P(k)*(EXPFU .03A16A75/ALMIK) )*LDGF(TMR!K)/TM) ) ) 

10 DNS I TY * PRESS/I2. 8704«TM) 

GO TO 13 

C 

C CONTROL COMES HERE FOR ISOTHERMAL LAYERS 

11 TM * TMRIK) 

PRESS” REF P(K)»EXPF(-0.034i6475*H I NC/TMRI K ) ) 

GO TO 10 
C 

C CONTROL COMES HERE FOR EXTREME ALTITUDES 

12 PRESS * 0.0 
DNS I TY * 0.0 
TM * 2000. 

13 RETURN 
C 

C END OF THE FORTRAN STATEMENTS. 


AL 

A2 

A3 

A4 


REM 

REM 

REM 

REM 

REM 

REM 

REM 

REM 

REM 

REM 

ORG 

REL 

DEC 

DEC 

DEC 

DEC 

DEC 

DEC 

DEC 

DEC 

REM 

END 


THIS IS THE SAP PROGRAM WHICH LOADS ICAO OATA INTO MACHINE. 

THE 170 IN ORG 170 WAS FOUND BY SUBTRACTING 10 FROM THE DEC LOCATION 
□F REF P (FROM SAP LISTING UF ICAO, THIS WAS FOUND TO BE 180). 

THUS, 180-10*170. 

Al IS REF Pill) 

A2 IS ALM(ii) 

A3 IS TMRUl ) 

AA IS TABLE Hill) 

170 


1.01E-8, l.A77E-8,6» 19E-7, 1. 451E-5 , l .815E-3,2.452E-2» 5. 832E-1 

1.20AA,2A.B86, 226. 32, 1013.25 

0.0,0.0005,0.0058,0.01 ,0.035,0.0,-0.0039,0.0,0.003,0.0 
-0.0065 

0.0, 15 37.86,812.86,322.86, 196.86, 196.86,282.66,282.66,216.66 

216 . 66 , 288.16 

3000000.0. 300000.0.175000.0. 126000.0.90000.0.75000.0.53000.0 

47000.0. 25000.0.11000.0.0.0 
END OF THE SAP STATEMENTS. 
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SUBROUTINE THRUST 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c 

c 


THIS ROUTINE COMPUTES X,Y,Z THRUST ACCELERATIONS. THE THRUST VECTOR IS 
ASSUMED COINCIDENT WITH THE LONGl TUNDI NA L AXIS OF THE VEHICLE, WHICH IS 
ORIENTED TO THE RELATIVE WIND VELOCITY BY THE ANGLE OF ATTACK [ALPHA) AND 
THE ROLL ANGLE (BETA). ALPHA IS ASSUMED TO BE A QUADRATIC FUNCTION OF TIME 
WHEREAS BETA IS ASSUMED TO BE CONSTANT. 

REVOLV IS THE EARTHS ROTATION RATE IN RADIANS/SEC ( 7. 2921 1585E-5 ) AND THE 
FACTOR 8589934592.= 2**33 IS REMOVED TO PREVENT OVERFLOW. 


COMMON C 


DIMENSION FORCE 

(3), P AR ( 3 ) , C ( 30 ) , VA TM ( 3 ) , P ( 

3) , AQ( 5) , INDI 3) 

EQUIVALENCE 

SI HP i C( 5)), 

F LOW , C ( 33) 

),( FORCE.CI 525) ) 

2 ( 

RMASS, C( 131 ) ) , 

PAR,C(798) ) , 

RSQRD,C( 567) 

) , (COSBET, Cl 599) ) 

3 { 

VX , C ( 472 ) ) , 

IND,C( 791) J , 

X,C< 200) 

> , ( S I NALF , C( 569)) 

41 

VY,C (473) 1 , 

TIME ■ C ( 1 38) ) , 

Y, C ( 20 1 ) 

) , (SINBET, C(568) ) 

5 ( 

VZ , C ( 474) ) , 

COSALF ,C( 575 ) ) , 

Z , C ( 202 ) 

) , (REVOLV, C(250) ) 

6 1 ALPHA , C ( 564 ) ) « 

I PMAGN ,C( 574) ) , 

P , C ( 5 7 1 ) 

>, ( R ATMOS , C ( 248 ) ) 

7 ( 

BETA , C 1 565) ) , 

[ VQSQRO ,C 1 481 ) ) , 

R , C ( 442 ) 

),( VATM , C ( 477 ) ) 

81 

VQ , C ( 480 ) ) , 

[ PUSH ,C < 34)) 




SINBET = SINF { BETA ) 

COSBET = C05F (BETA) 

VATM [ 1 1 * VX+REVOLV*Y 
VATM (2) =VY-R£VOLV*X 
VATM { 3 ) *VZ 

CALL CONVTUVATM, AQ) 

ALPHA = QUADITIME, l)/57. 29577951 
SINALF=SINF(ALPHA) 

CQSALF*COSF( ALPHA) 

DO I J 1* l , 3 
J2* I ND ( J 1 ) 

J 3= I NO I J 2 ) 

1 P(Jl> = (VATM[ J2)*AQ( J3 )-VATM( J3) *AQ( J2) )/ 0 5899 34592 . 

PMAGN* SQRTF(P(l)*P(i)*P(2)*P(2)*P(3)»P(3) J 

TDPMAG * PUSH/RMASS/PMAGN 
R4 * SINBET/VQ 

K 5 - COSALF/AQm 

DO 2 J 1 * l , 3 
J2 * I ND [ J 1 ) 

J3* 1 NO I J2) 

PARI J1 >*P{ J2) *VATM( J3) -P< J3 J *VATM( J2) 

2 FORCEIJl) * TDPMAG* (SINALF* (COSBET*P( Jl ) +R4*PAR( Jl ) ) -R5* { P ( J2 > *AQ 

I ( J3)-P( J3)*AQ( J2) ) ) 

RETURN 


END OF THE FORTRAN 


STATEMENTS. 


C 

C 

C 

C 

c 

c 

c 

C 

c 

c 

c 

c 


c 


c 

c 


c 


SUBROUTINE AERO 

SUBROUTINE AERO COMPUTES THE LIFT AND DRAG ACCELERATIONS. AS IN SUBROUT- 
INE THRUST, THESE VECTORS ARE REFERENCED TO THE RELATIVE WIND VELOCITY. 
COEFFICIENTS OF LIFT, INDUCED DRAG, AND DRAG AT ZERO ANGLE OF ATTACK ARE 
ASSUMED TO BE FUNCTIONS OF MACH NUMBER AND ANGLE OF ATTACK. TABLES OF 
CD I / CL * *2 , CL/SINI ALPHA), AND CDO ARE ASSUMED AS FITTED QUADRATIC EQUAT- 
IONS IN THE COEFN ARRAY. GASFAC IS THE SQRTF ( SPEC IF IC HEAT RATIO * STAND- 
ARD ACCELERATION OF GRAVITY * UNIVERSAL GAS CONSTANT). FOR EARTH, GASFAC* 
20. 064881 (METERS / SEC / KELVIN DEGREE). 


COMMON C 

DIMENSION C ID , VATM (3) , P(3> , X I F T ( 3 ) , DR AG( 3 ) , PAR ( 3) 


EQUIVALENCE 
1( Q VAL , C ( 794 ) 
Z( BETA , C ( 565 ) 
3< PHIP , C ( 462 ) 
4 ( XIFT ,C{528) 
5 ( DRAG ,C ( 53 l ) 
6IC0SALF »C ( 575 ) 
7 ( PAR , C ( 798 ) 


,( AREA, Cl 35) 
,( P MAGN , C ( 5 74 ) 
, (VQSQRD,C(481) 
,( RMASS ,C ( 131) 
,( VATM ,C 1477) 
,( VMACH,CI47l) 
,( CD , C l 7 97 ) 


, ( 

TIME, C( 138) ] 


,( 

TM, C ( 467) ] 


, ( 

VU, C( 480} ! 


, ( 

R, C ( 442 ) 1 


, ( 

P,C< 571 ) ] 


, < 

ALPHA, C(564) ] 


, ( 

CL, Cl 796) ] 



( ON S I TY , C ( 460 ) ) 
( S ! NALF ,C l 569 ) ) 
(SINBET, C(568>) 
< CD1 , C ( 795 ) ) 
(GASFAC, C( 458) ) 
(COSBET ,C(599) ) 


QVAL-O. 5*DNSI TY*VQSQRO* ARE A/RMASS 
VMACH=SQRTF( VQSQRD/TM)/ GASFAC 


COMPUTE THE X,Y,Z COMPONENTS OF LIFT. 

IF (ALPHA) 2,1,2 

1 CL * 0.0 
CD I *0. 0 
XIFT(l) * 0. 

XIFT 12) = 0. 

XIFT ( 3 ) = 0. 

GO TO 4 

2 CL * QUAD(VMACH,2)*SINALF 
DO 3 K= 1 ,3 

3 XIFT(K) * Q VAL*CL/PMAGN* (SI NBET*P AR ( K) /VQ+COSBET *P ( K ) ) 
CDI=QUAD(VMACH,3)*CL*CL 
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C COMPUTE THE X,Y ,2 COMPONENTS OF DRAG* 

4 CD « CD I +QUAD ( V MACH *4 ) 

DO 5 K = 1 » 3 

5 DRAG ( K ) * -CD*QVAL*VATM(K)/VQ 
RETURN 

C 

C END OF THE FORTRAN STATEMENTS, 


SUBROUTINE OBLATE 
C 

C THIS SUBROUTINE COMPUTES THE OBLATENESS ACCELERATIONS (OBLAT) DUE TO AN 
C AXIALLY SYMMETRIC EARTH. THE 2ND AND 4TH SPHERICAL HARMONIC COEFF, ARE 
C OBLATJ AND OBLATK, RESPECTIVELY. OBLATJ, OBLATK, RESQRD, AND THE CONSTANTS 
C CON ARE LOADED BY STDATA. 

C 

COMMON C 
C 

DIMENSION RB(3), OBLATI3), C0N(9) 

C 

EQUIVALENCE 

II CONfC I 576) ) , ( R tC ( 442) 1 $ I GK2M, C ( 469) ) , l RSQRD,C ( 567) ) , 

2 { RESQRD t C ( 40 >>,( OBLAT J , C ( 38) ) f (OBLATK, C( 39)), ( RB,C(200)) f 

3 ( OBLAT, C(534) ) 

C 

Z2DVR2 X RB(3)*RB(3)/RSQRD 
REDVR=RESQRD/RSQRD 
DO l K * 1 » 3 

l OBLAT(K)=RB(K) »RED VR*GK2M*5 . O/R/RSQRD* ( OBL ATJ* ( Z 2DVR2-CON I K ) >♦ 

1 OBLATK#REDVR*( Z20VR2* ( CQN( K+3 )-2. I *Z2DVR2 ) -CON I K+6 ) ) ) 

RETURN 

C 

C END OF THE FORTRAN STATEMENTS. •*#**#*< 


SUBROUTINE EPHMRS 
C 

C SUBROUTINE EPHMRS IS CALLED TO COMPUTE THE POSITIONS OF THE PERTURBING 
C BODIES RELATIVE TO THE VEHICLE AND, FROM THESE, THEIR PERTURBING ACCELERA- 
C TIONS UPON THE VEHICLE. OCCASIONALLY THIS ROUTINE IS CALLED FOR THE PURPOSE 

C OF TRANSLATING THE ORIGIN IN WHICH CASE { TRSF ER= 1 ) THE RELATIVE VELOCITIES 

C ARE ALSO CALCULATED. IF A BODYS POSITION IS TO BE COMPUTED FROM AN ELLIPTIC 
C APPROXIMATION SUBROUTINE ELIPSE IS CALLED. OTHERWISE, THE POSITION WILL BE 
C CALCULATED IN EPHMRS FROM THE PRECISION TAPE EPHEMER IS. THE 00 19 LOUP 

C ENCOMPASSES ALMOST THE ENTIRE EPHMRS SUBROUTINE AND , IN EFFECT, ELIPSE TOO. 

C 

COMMON C 
C 

DIMENSION QX ( 3 ) , I BODY ( 8) »EFMRS( 7) t XP( 3,8) » K B ( 3 * 8 ) , RREL ( 8) ,NEFMRS 
l (0) »TDATA( 18,7) »CF ( 6» 3 » 7 ) , T I M I 7 J , TDEL ( 7 ) , BMASS < 8) , XDOT ( 3, 8 ) ,C ( 1) 

C 

EQUIVALENCE (QX ,C(522)),l I BODY , C ( 425) ) , ( MBODYS , C ( 441 ) ) , 

1 { EFMRS ,C(4I0}}»(XP ,C(i76)),{RB , C ( 200) ),( RREL ,C(442>), 

2 ( NEF MRS » C ( 433 ) ) , ( TRSFER ,C(224) ) , ( TABL T,C(252) ) , ( DTOFF J,C< 31 )) , 

3 ( TDAT A , C ( 276 ) ) , (CF ,C ( 276 ) ) , ( T I M , C ( 585 )),( TDEL ,C(592)), 

4 ( BMASS »C ( 4 17 H , ( SQRDK ,C ( 468 ) ) , ( XDQT , C ( 498 ) ) , ( L ENGTH, C ( 2 57 ) ) , 

5 ( AU,C ( 461 ) 1 • ( IBFfFIB) 

C 

C PART 2. SET INDEXS, FIND POSITION IF ELLIPSE IS USED (NEF MRS * 20 OR UP). 

DO 19 JB=i, MBODYS 
JB1 - JB-H 
IBF * IBODYUBl) 

IB = XABSF(IBF) 

IF (NEFMRSt J6)-20J 2,2,1 

1 CALL ELIPSE (JBU 

IF (TRSFER) 12,12,17 
C 

C PART 3. TAPE EPHEMER I S IS TO BE USED. FIND DIFFERENCE (OT) BETWEEN 

C CURRENT PROBLEM TIME ( DTOFF J + T ABL I ) AND MIDPOINT TIME (TIM) OF CURRENTLY 

C STORED TAPE DATA. THEN SEE IF CURRENT DATA IS OKAY. TDEL = TIME INTERVAL 

C ON EITHER SIDE OF TIM FOR WHICH CURRENT DATA IS GOOD. 

2 DT - TABL T - (TIMIJB) -DTOFF J I 
IF (ABSF(DT)-TDEL(JB) ) 10,10,3 

C 

C PART 4A. CURRENT DATA NOT OKAY. READ IN NEXT DATA SET. IF DT IS 

C BACK UP THE TAPE 2 RECORDS BEFORE READING. 

3 IF (DT) 4,5,5 

4 BACKSPACE 3 
BACKSPACE 3 

5 READ TAPE 3, (C(J), J*8051,8071» 

LYE = 8051 

C 
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C PART 4b - IF THIS DAT* IS FOR A BODY IN THE BNAME LIST, STORE IT. 

C ( IF NOT STORED, WE MIGHT HAVE TO RETURN FOR IT,) IF ELLIPSE DATA IS 

C PROVIDED FOR THE BODY FOUND, BY-PASS THE TAPE DATA AND READ IN NEXT SET. 

DO 7 J = 1 , MBODYS 
S CLA C(LYE) 

S CAS EFMRS(J) 

S TRA *7 

S TRA *6 

S TRA *7 

6 IF l NEFMRSI J)-2Q ) B,8,3 

7 CONTINUE 
GO TO 3 

C 

c PART 4C , MOVE THE DATA INTO PLACE AND THEN GO BACK AND SEE IF IT IS OKAY. 

8 T I M ( J 1 » CtLYE+l) 

TDELU) = C ( L YE + 2) 

DO 9 J J- 1 , 1 8 

TDATA ( J J , J ) = C ( J J+8053 ) 

9 CONTINUE 
GO TO 2 

C 

C 

C PART 5. CURRENT DATA IS OKAY. GET POSITION FROM THE POLONQMIAL 

C P = A ♦ BX + CX»*2 ♦ DX**3 + EX**4 + FX**5. 

10 DO 11 K= 1 , 3 

XP(K, JB1 J = CF( 1,K,JB) 

DO II KT = 2 , 6 

XP ( K, JB 1 ) * XP(K, JB1 ) * DT +CF(KT,K,JB) 

11 CONTINUE 

IF ( TRSFER ) 12,12,15 
C 

C PART 6. COMPUTE DISTANCE FROM REFERENCE AND FROM ROCKET . 

12 DO 13 K= 1 , 3 

XPtK p JB1 ) = XPIK.IBJ +XP(K»JB1 1 *S IGNF ( AU , F I B ) 

13 RBIKi JBl I* RB t K , 1 ) - XP(K,J01) 

C 

C PART 7. COMPUTE PERT JRBI NG ACCELERATIONS (QX). 4194304=2**22 IS REMOVED 

C TO PREVENT OVERFLOW. 2048=2**11 AND 8509934592=2**33 RESTORE THE SCALE. 

PRSQRD » (RB( 1 , JB1 ) **2 * RB(2,JB1)**2 ♦ RB ( 3, JB 1 ) **2 1 /4 1 94304 . 

RRELL = 5QRTF ( PRSQRD ) 

RSQRD = t XP (1 , JB 1) **2 ♦ XP(2,JB1)**2 ♦ XP 1 3 , J B 1 ) **2 ) /4 1 94304 . 

RCUBE = RSQRD * SQRTF ( R SQRD1 
PRCUBE * PRSQRD * RRELL 
RREUJB1) = RRELL* 2040. 

DO 14 K= 1 , 3 

14 QX(K)=SQRDK * BMASSUB1) * {( XP( K, JB1 ) /RCUBE ) ♦ RB I K , JB 1 ) /PRCUBE ) / 

1 8589934592. ♦ QX(K) 

GO TO 19 
C 

c PART 8. COMPUTE VELOCITY FROM V = B + 2CX ♦ 3DX»*2 + 4EX**3 ♦ 5FX**4 

C AND FROM REFERENCE BODY VELOCITY ( XDOTC IB ) ) . 

15 DO 16 K= 1 , 3 
XDOTlK.JBl) = 0. 

DO 16 KT= 1 , 5 

16 XDOT ( K , JB 1) = ( XDOT ( K , JB 1 ) * DT + CF(KT,K,JB) *FLOATF l -KT*-6 ) ) 

17 DO IB K= 1 , 3 

18 XDOT ( K , JB 1 ) = XDOT ( K, IB) + X DOT ( K, JB 1) *S I GNF ( AU/86400. 0 , F I B ) 

GO TO 12 

19 CONTINUE 

CALL DUMP (4,C, LENGTH) 

RETURN 

C 

C END OF THE FORTRAN STATEMENTS. *•*.***. 


SUBROUTINE ELIPSE ( JBU 
C 

C THIS SUBROUTINE IS CALLED FROM EPHMRS TO COMPUTE THE POSITION OF A BODY 

C USING APPROXIMATE ELLIPTIC DATA. THE VELOCITY IS ALSO COMPUTED IF THE 

C ORIGIN IS BEING TRANSLATED ( TRSFER= 1 . 0 ) . THE ELLIPSE DATA IS READ FROM 

C INPUT CARDS AND ORGANIZED IN SUBROUTINE ORDER. TPD IS TIME SINCE PERIHELION 

C PASSAGE, ZM IS MEAN ANOMALY, U IS ECCENTRIC ANOMALY, E IS ECCENTRICITY. 

C 

COMMON C 
C 

DIMENSION 

1 XP (3,8), 

2 E < 1), 

3 SOMEGA ( 1 > , 

4 PERIOD (1), 

5 COMEGA ( 1 ) 

C 


XDOT (3,8), 
5INCL m, 
PPJD (1), 
CINCL (Hi 


P (1) , 
SNODE (1), 
PPFRAC (U* 
C NODE (1), 
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EQU I VALENCE 

i( XDOT ,C (498) ) * (DTOFFJ , C( 3 1 ) ) , ( COME GA, C ( 28A ) ) , { CNODE , C \ 285 } ) , 
2( P , C ( 276 ) ) , ( E.CI277) ) , ( SOMCGA , C < 2 > , I SNODE , C ( 2 79 H * 

3 ( SINCL , C (280) ) , ( PPJD,CI 231 ) J , ( PP Fk AC , C < 2 62 ) ) , I P ER I DO » C l 2 63 ) ) , 
A{ CONSU i C ( 36 ) ) * ( TABLT ,C ( 252 ) ) « ( XP , C I 1 76 >) , I TRSFER , C < 22A > ) , 
51 C I NCL * C ( 286 ) ) 

K * 18* C JB 1-2 ) ♦ l 

TPD = {OrOFFJ-PPJD(Kn + (TABLT-PPFRAC(K} ) 

IN = 6.26318533/PERIOD< K) 

ZM = ZN*M0DF ( TPD, PER ICO I K } ) 

GET THE SINE(SINTRU> AND THE COSINE (COSTKUJ OF THE TRUE ANOMALY 
BY ITERATING KEPLERS EQUATION. THEN COMPUTE X,Y,Z (XP>. 

U = ZM + EIK)*SINF(ZM)+0.5*E(K)**2*SINH2.0*ZM) 

DO 1 J - 1 » 10 

DELM = ZM-U+E(K)*SINF(U) 

DELU = DELM/I 1 . -E ( K ) *COSF ( U ) ) 

U = U+DELU 

IF ( ABSF (DELM)-CONSU) 2,2il 

1 CONTINUE 

2 COSU = COSF(U) 

DENGM = 1 • -E I K ) *COSU 
COSTRU = I COSU-E(K) I/DENOM 
R = P { K ) / { 1 • +E { K ) *COS TRU ) 

SI NTRU= SuRTF I 1 . -E ( K ) **2 ) *S I NF ( U ) / DENOM 
SINV = S 1 NT RU*C OMEGA ( K ) + COSTRU* SOMEGAI K ) 

CUSV = CUSTRU*COME GA ( K) -S I NTRU* SOMEGA ( K ) 

XPU, JBII = R* ( COS V * C NODE! K)-SINV*S NODE < K)*CINCL (K) ) 

XP(2,Jbi) = R*(COSV*SNODEiK)+$lNV*CNODE(K>*ClNCL(K) ) 

XP(3»Jbl) * R*$INV*SINCLIK) 

IF (TftSFER) 3, A, 3 

COMPUTE THE VELOCITIES FOR TRANSFER OF ORIGIN. 

3 EX = t ( K ) *SOMEGA(K) +SINV 
FX = EU)»COMEGA(K)+COSV 

CFACT = ZN*P(K) /( SQRTFI { i . 0-E ( K > * *2 ) * * 3 > > 

AX = EX*CNODE<K)*FX*SNODE(K)*CINCL(K J 
BX = FX*CN0D£<K|*C1NCL(K)-EX*SN00E(K> 

XDOT ( 1 1 JB L ) * -AX*CFACT 
XDOT { 2 , J 6 1 ) * BX *C F AC T 
XDOT ( 3 • JB L ) = FX*CF ACT*SINCL (K I 
A RETURN 

END OF THE FORTRAN STATEMENTS. 


SUBROUTINE CONVTUVX,A) 


THIS ROUTINE COMPUTES 


(1) ANGULAR MOMENTUM, A(A) 

(2} ANGULAR MOMENTUM SQUARED, A{5) 

{3) X , Y » Z COMPONENTS OF ANG. MOM., A ( 1 J , A { 2 ) , A ( 3 ) 
(A) VELOCITY, VX(A> 

(5) VELOCITY SQUARED, VX(5) 


COMMON C 
C 

DIMENSION A(5) , VX( 5) ,X( 3) , IND( 3 J 
C 

EQUI VALENCE (X,CI200)),(IND,C1791}) 

C 

DO 1 J 1=1 , 3 
J2= I ND l J 1 ) 

J 3= I ND ( J2 ) 

1 A(J3) = X{J1)*VX(J2)-XU2)*VXU1J 
A( 5)=A(1)*AC i>*A(2)*A(2)+AI 3)*A(3) 

A { A J = SQRTF I A ( 5 ) ) 

VX(5i=VX< 1)*VX( IUVXI2) *VX(2H-VX(3)*VX(3) 
VX(A} = SQRTF I V X I 5) ) 

RETURN 

C 

C END OF THE FORTRAN STATEMENTS. 
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SUBROUTINE C0NVT2 


THIS ROUTINE CONVERTS RECTANGULAR COORDINATES INTO ORBIT ELEMENTS. 
RECTANGULAR COORDINATES- POSITION COMPONENTS , X , AND VELOCITY COMPONENTS , VX. 
ORBIT ELEMENTS - tl) ECCENTR I C I T Y , E (4) INCLINATION, ZINCL 

(2) ARG. OE PERICENTLR, OMEGA (5) MEAN AMOMAL Y , ZM A 

( 3 ) LONG. OF ASCENDING NODE , ZNODE (6) StMILATUS RECTUM, P 


COMMON C 


DIMENSION CIl) ,VX(3) ,X(3) 


EUUI VALENCE 


1( 

A2 , C ( 5 59 ) ) , 

( OMEGA, C ( 228 ) ) 

2 ( 

A 3 , C ( 560) ) , 

( ZNODE ,C(229) ) 

3 { 

A l , C ( 56 l ) i , 

( ZINCL, C(230) ) 

4 [ 

P,C(232) ) , 

{ ZMA , C ( 2 3 1 ) ) 

5 ( 

R , C ( 442 ) ) , 

I S I NTRU ,C < 492 ) 1 

6 ( 

E , C ( 22 7 ) ) , 

( A , C ( 562 ) ) 


,( 

ASQRD.U 563) ) , 

( VX»C(472) ) • 

,( 

V,C( A75) 1 , 

( GK 2M , C ( 469 ) ) , 


VSORD.Ct A/6) ) , 

( EPAR , C ( 245 1 ) , 

• t 

TRUi C ( A03 ) ) , 

<TRSFER,C( 224) ) , 

, (C0STRU,C(493) ) , 

< X , C ( 200 ) ) , 


1 

2 

3 


4 

5 

6 
7 


P= ASQRD/GK2 M 

R = SQRTF ( X ( 1 ) * *2+ X < 2 ) **2 + X ( 33 »»2 I 

TRU=ARCTAN(A/GK2M*(X( I ) *VX (1) +X ( 2 ) * VX ( 2 ) +X { 3 ) *VX ( 3 )) , P-R ) 

IF { A2 ) 2,1,2 
ZNODE = 0.0 
GO TO 3 

ZNODE = ARC TAN { A2 , -A3 I 

ZINCL = ARCTANt SQRTFI A2**2+A3*«2) , Al) 

SNODE = S I NF ( ZNODE ) 

CNODE = CDS F ( Z NODE ) 

XT WOO = X ( l ) *CNODE *X ( 2 ) * SNODE 

YTWOD = X(3)*SINF (ZINCL) + COSF(ZINCL) * ( X ( 2 ) *CNODE- X ( 1 ) * SNODE ) 
OMEGA=ARCTAN( YT WUD , XTWOD ) -TRU 
E = SQRTF (ABSF( 1.4-p# { VSURD/GK2M-2 • / R ) ) ) 

EPONE = SORTF ( 1 .+E ) 

E2M1 = l.-E*C 

EPAR = SQRTF ( ABSF ( E2MI ) ) 

S1NTRU-SINF (TRU) 

CO STRU = COSF ( TRU) 

EPAS = SaRTFtABSFI l.-h) )*SINTRU/I UOfCOSTRU) 

ETHETA = E*SINTRU/( l . 0*E * COS TRU ) * EPAR 
IF { E2M i ) 5,6,6 

ZMA * LOGF ( ( E PONE + EPAS) /(EPONE -EPAS) ) - E THETA 
GO TO 7 

ZMA = 2.0*ARCTAN(tPA5, EPONE ) - ETHETA 
RETURN 


END OF THE FORTRAN STATEMENTS. 


FUNCTION ARCTAN (Y,X) 

THE FORTRAN II LIBRARY AT ANF { ♦ OR - Z= TAN ( THETA ) ) USES A SINGLE 
ARGUMENT WITH ITS SIGN TO GIVE THETA IN THE FIRST ( +Z ) OR FOURTH 
{ -Z ) QUADRANT. \ 

THE ARCTAN FUNCTION MAY BE USED IF ♦ OR - Z IS DERIVED FROM A 
FRACTION SO THAT ARCTAN (Y,X) = T AN- I { ( +0R- Y= S I N ( THETA > ) / ( +OR-X = 
cost THETA) )) . THUS THE ARCTAN (Y,X) GIVES THETA IN ITS PROPER 
QUADRANT FROM -180 DEGREES TO +180 DEGREES. 

IF (X) 2,1,2 

1 ARCTAN=S IGNF ( i . 57079632 *Y) 

GO TO 4 

2 ARCTAN=ATANF (Y/X) 

I F ( X ) 3 , i , 4 

3 ARCTAN=ARCTAN+S1GNF (3.1 4159265, Y) 

4 RETURN 

END OF THE FORTRAN STATEMENTS. 
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SUBROUTINE OUTPUT 


THIS IS THE ROUTINE WHICH FORMS THE BASIC DATA OUTPUT • BOTH ORBIT ELEM- 
ENTS AND RECTANGULAR COORDINATES ARE OUTPUT T ED* IF THE OBJECT IS NOT WITH 
IN AN ATMOSPHERE (PRESS=0.) f ONE LINE OF DATA IS DELETED. LIKEWISE, 

ONLY THOSE PERTURBING BODIES PRESENT HAVE THEIR DISTANCES OUTPUTTED. 


C 


C 


C 


c 


s 


c 

c 


c 

c 


c 

c 


COMMON C 


DIMENSION 




l 

RREL 

( 8) , 

ORBELS (6), 

C (1), 

2 

BNAME 

(8) , 

RB ( 3 , 8 ) , 

DIRC0S<3»8), 

3 

VAR 

(4) 




EQUIVALENCE 
I< TABLT ,C(252) > 

21 E » C ( 227 ) ) 

31 ZM A, Cl 231 ) ) 

4{ V«CI475M 

5( RREL » C < 442 ) ) 
6 I STEPGO t C ( 1 0 1 ) } 
7<DIRC0S,C(176M 
8 ( MBODYS , C ( 44 1 ) ) 
9 ( S I N TRU , C ( 492 ) ) 
EQUIVALENCE 
II ALT , C ( 463 ) I 

21 VQ,C(480)) 


,< TIME • C C 138) ) 
,< OMEGA, CI220) ) 
,< P,CI 232)1 
,( VX,C(472)I 
,( X , C ( 200 ) ) 
,1 DELT , C I 2561 ) 
, (0R8ELS,C(227) ) 
, ( NBODYS , C { 489 ) ) 
, (CQSTRU,C(493> ) 

,( VATM1 ,C(477) ) 
,( PSI,C(462J) 


I ST EPNO, C ( 102 ) ) 
I ZNODE , C I 229 ) ) 
I RB, C ( 200) ) 

I V Y » C ( 473) ) 
( Y,C(20D) 

I RMASS, C ( l 3 1 ) I 
( IMODE, Cl 28)) 
I DT OFF J , C I 31 ) ) 
I RE VS , C ( 490 ) ) 

( V ATM2 , C { 478 ) ) 


IBNAME , C ( 402 ) ) 
I Z INCL * C ( 2 30 ) ) 
I TRU,C(483)J 
( VZ , C ( 474 ) ) 
I Z , C I 202 ) ) 
I ALPHA, C(564) ) 
I PRESS, CI466J ) 
I A , C ( 562 ) ) 
(LENGTH, CI257) ) 

f VATM3 , C ( 479 ) ) 


PATHANF I VX, VY,VZ) - ATANFI ( X*VX+Y*VY*Z*VZ)/A)*57.2957795l 

DAYJ * IDT OF FJ-2.4E6) ♦TABLT 

ALPHA1 = ALPHA*57. 29577951 

REV = REVS ♦ ARCTANI-Y,-X)/6. 28318532 ♦ .5 

CALL C0NVT1IVX,C(559) ) 

LXD IMODE, (IMODE) 

GO TD (2,1,1) , IMODE 

1 C0DE=6HRECTAN 
CALL CON VT 2 
GO TO 4 

2 DO 3 K = i i 6 

3 ORBELS(K) = CIK+13L) 

CODE -5H0RB I T 

TRU=ARCTAN( SINTRU, CUSTRU) 

4 PSI = PATHANF (VX,VY,VZ) 

WRITE OUTPUT TAPE 6, 1 1 , STEPGO , STEPNO, E, OMEGA, V , RREL (1) , BNAME ( 1 ) , 
IC0D6 , 1 MODE, TIME, P, TRU, VX,X, RMASS, DAY J,ZMA, ZNODE, VY,Y, REV, ALPHA l, 
2PSI »ZINCL»VZ,Z,DELT 


IF WITHIN AN ATMOSPHERE COMPUTE DRAG, LIFT, G, ETC., AND PRINT EXTRA LINE. 
IF (PRESS) 5,7,5 

5 J = 0 

DO 6 1-1,4 
J ' J + 3 

6 VAR ( I ) = SQRTF(ClJ+525)**2+C{J+526)»»2*C(J+527)**2)*RMASS/9.80665 
G = VAR ( 4 ) /RMASS 

CALL C0NVT1I VATMi,C(559) ) 

PSI - PATHANF (VATMlfVAT M2, VATM3) 

WRITE OUTPUT TAPE 6 , 1 2 , ALT , PSI , VAR ( 2 ) , VQ ,G, VAR< 1 ) 

IF PERTURBATING BODIES ARE PRESENT, FIND THEIR DISTANCES AND PRINT THEM. 

7 IF(MBODYS) 8,10,8 

8 00 9 J = 2 , NBODYS 
DO 9 K=1 ,3 

9 DIRCOS(K,J) = -RBI K, J I/RREL I J) 

WRITE OUTPUT TAPE 6,13, 

1 ( 8 NAME ( J ) ,RREL( J) tDIRCOSll, J) , DIRC0S(2,J),DIRC0S(3,J),J=2, NBODYS) 

10 CALL DUMP(2,C, LENGTH) 

RETURN 

11 FORMAT ( 6H05TEP-F5. ,2H +F4.,4X,13H ECCENTRIC IT Y=1PG15.8,7H OMEGA-G l 5 
1 . 8 i 4H V=G15. 8 , 3H R=G15.8,7H REFER- A6, IX ,A6 , I2/6H T IME= IPG14. 7,14 
2H SEMILATUS R.=G15.8,7H TRU A-G15.8,4H VX-G15.8,3H X-G15.8,7H RMAS 
3S-G1 5. 8/9H JDAY- 240PF 1 0. 4 , 1 5H MEAN ANOMALY- 1PG1 5. 8 , 7H N0DE-G15. 

48, 4H VY-G15.8,3H Y=G15.8,7H REVS.-G 15. 8/6H ALFA-G14. 7, 14H PATH A 
5NGL E =G l 5 . 8 , 7H I NC L=G 1 5 . 8 , 4H VZ=G15.8,3H Z-G15.8,7H DELT-G15.8) 

12 FORMAT ( 6H AL T . = IPG l 4. 7 , 14H R PATH ANGLE-G1 5 . 8, 7H DRAG-G15. 8, 4H VR 
l = G 15 • 8 » 3H G = G15.B,7H LIFT*G15.8I 

13 F0RMAT(2(1X,A6,3H R- i PG14. 7 , OP 3F 10. 6, l IX ) > 


END OF THE FORTRAN STATEMENTS. 
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SUBROUTINE DUMP ( IDENT » DATA , LENGTH ) 


C 

c THIS subroutine will dump in g type format a variable number of consecutive 

C WOROS, BEGINNING AT A SPECIFIED LOCATION. DUMP OCCURS WHEN THE FOLLOWING 

C CONDITIONS ARE SATISFIED 
C 

C A) IDENTIFICATION NUMBER (IDENT) = AN INPUT DUMP NUMBER ( NDUMP ) . 

C B) DUMP NUMBER IDENT HAS BEEN SKIPED NSKIP TIMES. 

C C) TOTAL NUMBER (TEST) OF DESIRED DUMPS HAS NOT BEEN EXCEEDED. (IF TEST 

C IS NEGATIVE, DUMP ALWAYS OCCURS). 

C 

C NOTE- IDENT = IDENTIFICATION NUMBER OF DUMP 
C DATA = STARTING LOCATION OF DUMP 

C LENGTH = NUMBER OF CONSECUTIVE WORDS TO BE DUMPED. (ZEROES COUNT BUT 

C ARE NOT DUMPEO ) 

C 

COMMON C 
C 

DIMENSION 

1 DATA (1), 16 (6), DATA6 (6), 

2 NSKIPN (A), NDUMP (A), NSKIP (A) 

C 

EQUIVALENCE 

L( TEST , C ( !)),( NDUMP, C(268) ), { NSK I P, C ( 2 72 ) ) 

C 

C PART i. TEST FOR OVERFLOW AND DIVIDE CHECK. 

IF DIVIDE CHECK 1,2 

1 ASSIGN 2 TO N 
WORDI = 6HDIVIDE 
W0RD2 = 6H CHECK 
GO TO 6 

2 IF ACCUMULATOR OVERFLOW 3, A 

3 ASSIGN A TO N 
WORDI = 6HACC OV 
WORD2 = 6HER FLO 
GO TO 6 

A IF QUOTIENT OVERFLOW 5,8 

5 ASSIGN 8 TO N 
WORDI = 6HMQ OVE 
WORD2 = 6HR FLOW 

6 WRITE OUTPUT TAPE 6 , 7 , WORDI , W0RD2 , I DENT 

7 FORMAT! 1H02A6, I8H IDENTIFICATIONS) 

GO TO N, (2, A, 8) 

C 

C PART 2. DETERMINE IF DUMP MAY OCCUR. 

8 IF (TEST) 15,26,9 

9 DO 12 1 = 1, A 

IF ( I DENT-NDUMP ( I ) ) 12,10,12 

10 IF (XABSF(N$KIP(I))-NSKIPN( I)) I3,L3,LL 

11 NSKIPNl I ) - NSKIPNI I IH 

12 CONTINUE 
GO TO 26 

13 NSK I PN { I ) = 0 

IF ( NSK IP ( I ) ) 1A,15,I5 
1A NSK I P ( I ) = 0 


C PART 3. DUMP OCCURS. OUMP NON-ZERO WORDS AND THEN REDUCE TEST BY l. 

15 WRITE OUTPUT TAPE 6 , 2 3 , TEST , IDENT , L t NOTH 
K2 = 6 

J = 0 

16 DO 21 K=i,6 

17 J s J*1 

IF (J-LENGTH) 18,18,19 

18 IF { DATA ( J ) ) 20,17,20 

19 K2=K-1 

IF ( K2 ) 22,25,22 

20 DATA6 ( K ) =DATA ( J ) 

21 I 6 ( K ) * J 

22 WRITE OUTPUT TAPE 6 , 2A , ( 16 ( K l ) , UATA6 (K 1) , K 1= 1 , K2 ) 

23 FORMAT ( 12H0DUMP, TEST=F6. I , 18H IDENTIFICATION IS,IA, 20H, NUMBER 
L OF WORDS IS, 15) 

2 A FORMAT ( IX , I A , 1PG1 5. 8 , 5 ( I 7 , IPG 1 5. 8 ) ) 

GO TO 16 

25 TEST = TEST-1. 

26 RETURN 
C 

C END OF THE FORTRAN STATEMENTS. *** 
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on ci o o ooooonoon 


FUNCTION QUAD IX, IC) 


THIS ROUTINE COMPUTES ANY VARIABLE* QUAD, AS A QUADRATIC FUNCTION OF X. 
QUAD * A ♦ BX * CXX. THERE MAY BE SEVERAL SETS OF COEFF I ENTS, EACH SET 
BELONGING TO A PARTICULAR REGION OF X. THE COEFN ARRAY IS ARRANGED AS — 

XI, A1 , Bi ,C I , X2, A2, B2 ,C2 * X3, A3, B3 * C3, X4, 

WHERE A I , B 1 * C 1 ARE THE COEFFI ENTS TO BE USED FOR X BETWEEN XI AND X2,ETC. 
AND XI IS LESS THAN X2 , X2 IS LESS THAN X3, X3 IS LESS THAN X4, ETC . 

IC IDENTIFIES WHICH DEPENDENT VARIABLE, QUAD, IS BEING SOUGHT. 

ICCIIC) DEFINE THE STARTING LOCATIONS IN THE COEFN ARRAY FOR VARIABLES X. 


COMMON C 

DIMENSION C(10> , COEFN I 190} , ICCC 51 

EQUIVALENCE i ICC, 01230)1, ( COEFN, CI601) > 

I = ICCUC) 

l IF (X-COEFN(n) 2,3,3 
21 s 1-4 
GO TO l 

3 IFU-COEFNI 1+4) ) 5,5,4 

4 I * 1+4 
GO TO 3 

5 QUAD = COEFN! I ♦ 1 ) ♦ X* ( COEFN { 1+2 ) +X*COEFN( 1*3)) 
ICC < IC ) = X 

RETURN 

END OF THE FORTRAN STATEMENTS. 



REM 

SUBROUTINE EXADD (A,B,C) 




REM 

THIS ROUTINE WILL ADO IN DOUBLE 

PRECISION A QUANTITY C 

TO THE DOUBLE 


REM 

PRECISION VARIABLE A+B WHERE A 

IS THE MOST SIGNIFICANT 

PART AND B IS 


REM 

THE LEAST SIGNIFICANT PART. 




ORG 

0 




PGM 





PZE 

END+1,0,0 




PZE 





BCD 

IEXAOD 




PZE 

EXADD 




ORG 

0 




REL 




Q1 

SYN 

32700 



Q2 

SYN 

32701 



TEMPI 

SYN 

32 702 



TEMP2 

SYN 

32703 




BCD 

iEXADD 



EXADO 

CLA 

1.4 




5TA 

TOPI 




STA 

T0P2 




CLA 

2,4 




STA 

BDTI 




STA 

BDT2 




CLA 

3,4 




STA 

ARGl 



TOPI 

CLA 

• « 



ARGl 

FAD 

#* 




STQ 

Ql 



B0T1 

FAD 

• • 




STQ 

Q2 




FAD 

Ql 




STQ 

Ql 




STO 

TEMPI 




CLA 

Ql 




FAD 

Q2 




STO 

TEMP2 




FAD 

TEMP2 




FAD 

TEMPI 




STQ 

Ql 




FSB 

TEMP2 



T0P2 

STO 

• • 




STQ 

Q2 




CLA 

Ql 




FAD 

Q2 



B0T2 

STO 

• • 



END 

TRA 

4,4 




REM 

END OF THE SAP STATEMENTS. 




END 






SUBROUTINE BKFILE(N) 

C THIS ROUTINE SIMPLY BACKSPACES TAPE N ONE FILE. 

C 

S 1 CAL *1 

S STP *Z 

M * 10-N 

S BST 10, (M) 

S BST 10, (M) 

S 2 BST 10, (MI 

S NOP 

5 RTB 10, (M) 

S CPY OUO 

S TRA * 3 

S TRA 

S TRA *3 

S 3 BST 10, (MJ 

A RETURN 
C 

C END OF THE FORTRAN STATEMENTS. 


REM SUBROUTINE PONG(N) 

REM THIS ROUTINE FINOS THE SEGMENT N ON TAPE ANO LOADS IT IN THE CORE. 
REM IF SEGMENT N IS ALREADY IN THE CORE, CONTROL IS SIMPLY SWITCHED TO 
REM THE BEGINNING OF SEGMENT N. 

ORG 

PGM 

PZE END* 1 , 0, 0 
PZE -1 
BCD 1SELPGM 
PZE SELPGM 
BCD 1PING 
PZE PING 
BCD 1P0NG 
PZE PONG 
REL 
ORG 

SPCO PZE DEC IS TOTAL RECORDS, ADDRS IS THIS RECORD, SET BY PING-PONG. 

PING TSX SPC A , 1 
SELPGM CLA 1,4 
PONG CLA 1 , A 
STA 

CLA *■ 

SSP 

TZE PING PROGRAM NUMBER TOO SMALL. 

SUB SPCO COMPARE WITH TOTAL 

TPL ERR PROGRAM NUMBER TOO LARGE. 

ADD SPCO 
ARS 18 

STO COMMON DESIREO NUMBER IN ADDRESS. 

CLA SPCO PRESENT POSITION IN ADDRESS. 

ADM SPC7 AOO ONE IN ADDRESS. 

ANA SPCB SAVE AOORESS ONLY. 

SUB COMMON 

PAX ,1 NUMBER OF RECOROS TO MOVE. 

TXL SPCA,1,0 PROPERLY POSITIONED IF ZERO. 

TXH SPC1,1,1 CORE LOAO OK IF ONE. 

PXO 

TRA SELPGM GO TO THE TRANSFER TO BEGINNING OF PROGRAM. 

SPCl TMI SPC2 ADVANCE TAPE. 

BST T BACKSPACE TAPE. 

TRA SPC3 
SPC2 RTB T 

SPC3 T I X SPCl, 1,1 KEEP MOVING TAPE. 

SPC4 RTB T 

CPY 0 RIGHT POSITION NOW SO LOAD IT. 

TRA SPC5 

RE W T EOF, NEXT IN SEQUENCE IS FIRST RECORD. 

TRA SPC4 FALSE EOR. 

SPC5 CAL 0 

ANA SPCB AOORESS OF FIRST WORD IS REQUIRED NUMBER. 

SUB COMMON DESIRED NUMBER. 

TXH SPC 7 ,1,1 BYPASS CHECK ON SELPGM ENTRY. 

TZE SPC7 PROPERLY FOUND. 

SPC6 HPR 1,5 IMPROPER POSITIONING OUE TO MACHINE ERROR. 

TRA SPC6 
SPC 7 CPY 1 

TRA 0 GO TO LOADER. 

SPCB PZE -l IS ONEA. 

ERR 1 BCD 10P0NG 

ERR2 BCD 1FAIL. 

ERtt WTO 6 THIS IS THE ERROR PRINTOUT ROUTINE. 

CPY ERR1 
CPY ERR2 
IDD 

RDR l CALL MONITOR. 

CPY 0 
CPY 1 
END TRA 0 

T EQU 2 

COMMON SYN -1 

REM END OF THE SAP STATEMENTS. ••••». 

END 
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TABLE II. - PROGRAM CONTROL PARAMETERS 


Control 
vari ables 

COMMON 

location 

Possible values 

Setting 

Description of use 

ASYMPT 

543 

0.0 or 1.0 

Internal 

Normally equal to 0.0, set equal to 1.0 in SUBROUTINE EQUATE when Kepler’s equation 
fails to converge for e > 1, and then used to control branching In MAIN 2 for 
IMODE - 3. 

ATKN 

26 

Any ALP coded body 
name 

Input 

Contains name of body which Is to have an atmosphere. Causes SUBROUTINE AERO to be 
called In SUBROUTINE EQUATE If object Is within that atmosphere. 

CLEAR 

19 

Any value 

Input 

If CLEAR - 0, SUBROUTINE STDATA is called from MAIN If If CLEAR / 0, 
SUBROUTINE STDATA Is bypassed. STDATA clears C{4) to C(1300). 

CONSTU 

18 

>0, ~10" B to ~10 -2 
radian 

STDATA: 10“ 6 

Input 

Controls branching in SUBROUTINE EQUATE, which determines flow accurate eccentric 
anomaly will be computed by Kepler's equation. 

CONSU 

36 

>0, ~10" 5 to -10 -2 
radi an 

STDATA: 10‘ 6 

Input 

Similar to CONSTU except that it is used in SUBROUTINE ELIFSE for perturbing bodies 
Instead of object. 

DELMAX 

23 

Any number of seconds 

Input 

If WODOUT - 2 or 3, output is given only at intervals of DELMAX. 

EHEF 

37 

Any number 

STDATA: 10 -6 

Input 

Desired error value. Error control predicts step size such that E2 ~ EREF. If 
5 REE <0, It will be treated as +EREF; however, error data will be recorded and 
printed , 

ERLIMT 

17 

Any plus number 

STDATA: 3xl0 -e 

Input 

Maximum error value that allows step in question to be passed as good step. 
If E2 > ERLIMT, step is recomputed with smaller step size. 

ETOL 

25 

Positive number of 
order 0.01 

STDATA: 0.01 

Input 

If eccentricity Falls in region 1 ± ETOL and Integration Is in orbit elements. 

Integration mode Is switched to temporary rectangular until eccentricity falla out- 
side this region. 

PILE 

249 

Any pluB Integer 

Internal 

Set equal to 10.0 in SUBROUTINE ORDER If tape data Ib used to determine positions, 
velocities, and attractions of perturbing bodies. Then read as file number of 
tape 3 in MAIN 1. See TFILE. 

ICC(5) 

238-242 

Any fixed -point In- 
teger 

Input 

Internal 

Index of independent variable In COEJGl array used in FUNCTION QUAD. For each set of 
coefficients there is an ICC. Wey are aet at input time and are reset each time 
QUAD is called. 

IMODE 

28 

l,2,3,4,-l,-2,-3,-4 
(fixed point) 

STDATA: 1 

Input 

Internal 

Indicates Integration mode. Must agree with Input data (IT Input data Is rectangular, 
IMODE should equal 2 or -2). Values Indicate; 

1 - orbit elements -1 =* orbit elements, change to rectangular 

2 - rectangular variables -2 * rectangular, change to orbit elements 

3 * temporary rectangular -3 ■ orbit, change to temporary rectangular 

4 = Earth spherical change -4 ■ Earth spherical, change to orbit element 

to rectangular 

LENGTH 

257 

Any fixed-point in- 
teger 

Input 

Length of dump (i.e., number of words to be dumped). 

MODOUT 

103 

1,2,3, 4 
(fixed point) 

STDATA: 4 

Input 

Internal 

MODOUT - 1 Output every n th Step (n - STEPS) until TIME - WIN, then shift to mode 2 . 

m 2 Output at time intervals of DELMAX until TIME “ WAX. 

- 3 Output at time Intervals of DELMAX until TIKE - WIN, then shift to mode 4. 

- 4 Output every n tfi step until TIKE - WAX. 

HDUMP(4) 

268-271 

Any fixed-point In- 
teger 

Input 

If 1 In CALL DUMP (i, C, LENGW) command equals any number In NDUMP array, 
dump will be executed conditionally (see NSKIP). 

N3XIP(4) 

272-275 

Any fixed-point In- 
teger 

Input 

Causes skipping of NSKIP(i) dumps where N3KIP(i) corresponds to NDUMP (i ) . See SUB- 
ROUTINE DUMP. 

NPONG{5) 

11-15 

Any fixed-point in- 
teger 

STDATA: 2,1, , ,1 

Input 

NFONQ(i) refers to segment that is being called in statements CALL PONG (NPONG(l)). 
Control Is to beginning of segment. 

OBLATN 

27 

Any ALP coded body 
name 

Input 

If obl&teness effects are to be considered, loading a body name will cause SUBROUTINE 
OBLATE to be called from SUBROUTINE EQUATE when OBTAIN matches reference body. 

RECALL 

9 

Any value 

.. 

Input 

j 

If RECALL / 0.0, "starting" data will be restored from C(5) to C(115) in MAIN 1. See 
SAVE. 

SAVE 

8 

1.0, 2.0, or any 
other value 

Input 

! 

If SAVE * 1.0, "starting" data Trom c(5) to C(115) will be saved to be used later for 
another start requiring same data- If SAVE ■» 2.0, same thing happens, only before 
CALL INPUT (1) statement In MAIN 1. This saves result or previous Integration for 
future use. 

STEP GO 

101 

Any plus number 

Internal 

Total number of good steps. 

STEP NO 

102 

Any plus number 

Internal 

Total number of bad steps. Bad step does not pass error control teat. 

STEP MX 

20 

Any plus number 

STDATA: 100.0 

Input 

If (STEPGO + STEFNO) > STEPMX, problem terminates. 

STEPS 

21 

Any plus number 

STDATA: 1.0 

Input 

Used when MODOUT - 1 or 4. Output will occur at every n th step where n - STEPS. 

TAPE 3 

2 

0.0 or 3.0 

Internal 

Input 

If "working" ephemeris tape is to be made, TAPE 3 must be set equal to zero through 
input contained In SUBROUTINE TAPE. If no tape is to be made, or after tape is 
made, TAPE 3 la set to 3.0. 

TEST 

1 

Any integer 

Input 

Internal 

Total number of dumps. Initially set through input and thereafter decreased by one 
each time a dump occurs until TEST ™ 0. When TEST ■ 0.0 no more dumps will occur, 
ir negative value of TEST is loaded, there is no limit on number of dumps. 

TglTiR 

16 

Any plus integer 

[ STDATA: 1.0 

Input 

Selects which file of "working" ephemeris tape Is to be used. MAIN 1 positions 
tape in correct position by matching desired file number (TFILE) with code 
word (FILE) written at beginning of each file on tape. 

TMAX 

30 

Any number in seconds 

Input 

When TIME - TMAX control is switched to MAIN 1 to either read new input or end problem. 

WIN 

22 

Any number in seconds 

Input 

When TIME - WIN output mode Is changed. See MODOUT. 

TRSFER 

224 

0.0 or 1.0 

Internal 

Normally TRSFER = 0,0, but when origin is being translated TRSFER ■ 1.0 which 
causes SUBROUTINES EFHMRS and ELIPSE to compute velocities as well as positions. 

TTEST 

251 

Any number in seconds 

Internal 

When integration mode is changed to temporary rectangular, TTEST is set sb time at 
which program will begin checking for return to orbit elements. See MAIN 2, 
part 7D. 
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TABLE III. - BASIC OUTPUT FORMAT 


(a) Sample output 


STEP= 0. + 0. EC CENTRIC ITY= 1.00000000 0MEGA=-2. 64801353 

TIME= 0. SEMILATUS R.= 1 . 93844640E-09 TRU A= 3.14159262 

JDAY= 2437640.8350 MEAN ANOMALY= 0. NODE= 2.02516600 

ALFA= 0. PATH ANGLE= 89.9209976 INCL= 1.57079409 

ALT. =-0.1875000 R PATH ANGLE= 89.9209976 DRAG= 4. 99665982E-03 

SUN R= 1 . 4728028E 11 -0.261730 -0.885466 -0.383989 

V= 9. 99999976E-02 R= 6373346.50 REFER= EARTH RECTAN 2 

VX=-3 . 8 622435 9E-02 X=-2463371 . 37 RMASS= 150000.000 

VY= 7.90702742E-02 Y= 5043168.50 REV3.= 0.32231534 

VZ= 4. 74994606E-02 Z= 3019569.50 DELT= 6.00000000 

VR= 9. 99999976E-02 G= 1.49946962 LIFT= 0. 

MOON R- 3. 8293912E 08 -0.387660 -0.874846 -0.290456 


(b) Parameter identification 


FORTRAN code name 

Identification 

Output 

format 

mnemonic 

Internal 


STEP 

STEPGO, 

STEPNO 

Count of total number of successful integration 
steps to left of plus sign and count of fail- 
ures on right 

TIME 

TIME 

Time since beginning of Integration process, t, sec 

JDAY 

DAYJ 

Current Julian date 

ECCENTRICITY 

E 

Osculating orbit eccentricity, e 

SEMILATUS R. 

P 

Semllatus rectum of osculating orbit, p, m 

MEAN ANOMALY 

ZMA 

Mean anomaly of osculating orbit, M 

OMEGA 

OMEGA 

Argument of perl center, <x>, radians 

TRU A 

TRU 

True anomaly of osculating orbit, v, radians 

NODE 

ZNODE 

Equatorial longitude of ascending node of 
osculating orbit, n, radians 

INCL 

ZINCL 

Orbit Inclination referred to mean equator and 
equinox of 1950.0, I, radians 

ALFA 

ALPHA 

Angle between thrust and velocity, a, deg 

PATH ANGLE 

PSI 

Angle between path and local horizontal, deg 

V,VX,VY,VZ 

V,VX,VY,VZ 

Velocity and Its x,y,z components, V, m/sec 

R,X,Y,Z 

RREL(l), 
X, Y,Z 

Radius and Its x,y,z components, r, m 

REFER 

BNAME(l) 

Name of reference body, followed by integration 
mode, IMODE 

RMASS 

RMASS 

Vehicle mass, m, kg 

REVS. 

REV 

Revolutions past x-axis 

DELT 

DELT 

Step size for current step, h, sec 

ALT. 

ALT 

Altitude above oblate Earth, m 

R PATH ANGLE 

PSI 

Relative path angle, relative to Earth, deg 

DRAG 

VAR(2) 

Total drag force, D, kg 

VR 

VQ 

Velocity relative to rotating reference body 

G 

G 

Total Earth g’s acting on missile 

LIFT 

VAR(l) 

Total lift force, L, kg 

BNAME(i)R 

BNAME(i), 
DIR COS 

Vehicle to perturbing body distance, plus 

direction cosines 
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TABLE YI . - ASSUMED VALUES OF ASTRONOMICAL CONSTANTS 


Constant 

Assumed value 

FORTRAN 

name 

COMMON 

location 

Astronomical unit, m 

1 . 495X10 11 

AU 

461 

Gravitational constant, k 2 

1.32 452 139X10 20 

SQRDK 

468 

m 3 /(sec 2 )(sun mass units) 
Equatorial Earth radius 

4 .0680988 7 7X10 13 

RESQRD 

40 

squared, m 2 

Earth oblateness coefficient, J 

1.6238X10" 3 

OBLATJ 

38 

Earth oblateness coefficient, K 

6.4X10" 6 

OBLATK 

39 

Earth radii per AU 

4.26546512X10 -5 

ERTOAU 

a 3 

Day, sec 

86400 

SPD 

253 

Mass, reciprocal sun mass units: 
Sun 

1.0 

AMASS(l) 

881 

Mercury 

6,120,000 

AMASS(2) 

882 

Venus 

406,645 

AMASS(3) 

883 

Earth 

332,488 

AMASS(4) 

884 

Mars 

3,088,000 

AMASS(5) 

885 

Jupiter 

1047.39 

AMASS(6) 

886 

Saturn 

3500.0 

AMASS( 7 ) 

887 

Uranus 

22,869 

AMASS(8) 

888 

Neptune 

18,889 

AMASS( 9 ) 

889 

Pluto 

400,000 

AMASS (10) 

890 

Moon 

AMASS( 4 )/81 . 375 

AMASS ( 11 ) 

891 

Earth-moon 

AMASS( 4) + AMASS ( 11 ) 

AMASS ( 12 ) 

892 

Sphere-of-inf luence radii, m: 




Sun 

1.0X10 20 

RCRIT(l) 

911 

Mercury 

1.0X10 8 

RCRIT( 2 ) 

912 

Venus 

6.14X10 8 

RCRIT( 3 ) 

913 

Earth 

9.25X10 8 

RCRIT( 4 ) 

914 

Mars 

5.78X10 8 

RCRIT( 5 ) 

915 

Jupiter 

4.81X10 10 

RCRIT( 6 ) 

916 

Saturn 

5.46X10 10 

RCRIT( 7 ) 

917 

Uranus 

5.17X10 10 

RCRIT(8) 

918 

Neptune 

8.61X10 10 

RCRIT( 9 ) 

919 

Pluto 

3.81X10 10 

RCRIT( 10 ) 

920 

Moon 

1.60X10 8 

RCRIT(ll) 

921 


a Location relative to COMMON of subroutine TAPE (TAPE has a COMMON that is 
independent of all other subroutines ) . 
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TABLE VII. - LEWIS RESEARCH CENTER EPHEMERIS TAPE DATA 
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